/* ILapack.cs 
 * 
 * This provides our interface to LAPACK routines. 
 * Copyright (c) 2005, dnAnalytics. All rights reserved. 
 * 
 */

#region Using Directives
using dnAnalytics.Math;
#endregion

namespace dnAnalytics.LinearAlgebra.Native
{
    /// <summary>
    /// An enumeration of the available LAPACK providers.
    /// </summary>
    public enum LapackProvider
    {
        /// <summary>
        /// Provider for AMD's ACML.
        /// </summary>
        Acml,

        /// <summary>
        /// Provider for AMD's ACML compiled with Intel's FORTRAN compiler. A separate provider is needed
        /// because this version uses different export names than the other versions of ACML.
        /// </summary>
        AcmlIFort,

        /// <summary>
        /// Provider for the CLAPACK library.
        /// </summary>
        /// <remarks>ATLAS versions of the getrf, getri, getrs, potrf, potri, and potrs routines
        /// are used instead of the CLAPACK versions. ATLAS also provides the BLAS routines used by
        /// CLAPACK.</remarks>
        Clapack,

        /// <summary>
        /// Provider for INTEL's Math Kernel Library.
        /// </summary>
        Mkl
    }
    
    /// <summary>
    /// 
    /// </summary>
    internal interface ILapack
    {
        /// <summary>
        /// Computes the least squares solution to an over-determined system of linear equations, A*X=B, or the minimum norm
        /// solution of an under-determined system, where A is a general rectangular matrix of full rank, using a QR or LQ factorization
        /// of A or under determined linear system with full rank matrix.
        /// </summary>
        /// <param name="trans">If <b>trans</b> == 'N': the linear system involves <paramref name="a"/>. If <b>trans</b> == 'T': the linear system involves <paramref name="a"/>'.</param>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <b>m</b> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <b>n</b> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns of the matrix <paramref name="b"/>. <b>nrhs</b> &gt;=0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix. On exit, if <paramref name="m"/> >= <paramref name="n"/>, 
        /// <paramref name="a"/> is overwritten by details of its QR factorization as returned by Sgeqrf. Otherwise, <paramref name="a"/> is 
        /// overwritten by details of its LQ factorization as returned by Sgelqf.</param>
        /// <param name="lda">The leading dimension of <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b"> Dimension (<paramref name="ldb"/>,<paramref name="nrhs"/>). On entry, the matrix <b>b</b> of right hand side vectors, stored column wise. 
        /// b is <paramref name="m"/> by <paramref name="nrhs"/> if <paramref name="trans"/> == 'N', or <paramref name="n"/> by <paramref name="n"/> by <paramref name="nrhs"/>
        /// if <paramref name="trans"/> == 'T'. On exit, <b>b</b> is overwritten by the solution vectors, stored column wise: 
        /// <list>
        /// <item>If <paramref name="trans"/> == 'N' and <paramref name="m"/> &gt;= <paramref name="n"/>, rows 0 to <paramref name="n"/>-1 of <b>b</b> contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements <paramref name="n"/> to <paramref name="m"/>-1 in that column.</item>
        /// <item>If <paramref name="trans"/> == 'N' and <paramref name="m"/> &lt; <paramref name="n"/>, rows 0 to <paramref name="n"/>-1 of <b>b</b> contain the minimum norm solution vectors.</item>
        /// <item>If <paramref name="trans"/> == 'T' and <paramref name="m"/> &gt;= <paramref name="n"/>, rows 0 to <paramref name="m"/>-1 of <b>b</b> contain the minimum norm solution vectors.</item>
        /// <item>If <paramref name="trans"/> == 'T' and <paramref name="m"/> &lt;<paramref name="n"/>, rows 0 to <paramref name="m"/>-1 of <b>b</b> contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements <paramref name="m"/> to <paramref name="n"/>-1 in that column.</item>
        /// </list></param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <paramref name="ldb"/> &gt;= max(1,<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgels(char trans, int m, int n, int nrhs, double[] a, int lda, double[] b, int ldb);

        /// <summary>
        /// Computes the least squares solution to an over-determined system of linear equations, A*X=B, or the minimum norm
        /// solution of an under-determined system, where A is a general rectangular matrix of full rank, using a QR or LQ factorization
        /// of A or under determined linear system with full rank matrix.
        /// </summary>
        /// <param name="trans">If <b>trans</b> == 'N': the linear system involves <paramref name="a"/>. If <b>trans</b> == 'T': the linear system involves <paramref name="a"/>'.</param>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <b>m</b> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <b>n</b> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns of the matrix <paramref name="b"/>. <b>nrhs</b> &gt;=0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix. On exit, if <paramref name="m"/> >= <paramref name="n"/>, 
        /// <paramref name="a"/> is overwritten by details of its QR factorization as returned by Sgeqrf. Otherwise, <paramref name="a"/> is 
        /// overwritten by details of its LQ factorization as returned by Sgelqf.</param>
        /// <param name="lda">The leading dimension of <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b"> Dimension (<paramref name="ldb"/>,<paramref name="nrhs"/>). On entry, the matrix <b>b</b> of right hand side vectors, stored column wise. 
        /// b is <paramref name="m"/> by <paramref name="nrhs"/> if <paramref name="trans"/> == 'N', or <paramref name="n"/> by <paramref name="n"/> by <paramref name="nrhs"/>
        /// if <paramref name="trans"/> == 'T'. On exit, <b>b</b> is overwritten by the solution vectors, stored column wise: 
        /// <list>
        /// <item>If <paramref name="trans"/> == 'N' and <paramref name="m"/> &gt;= <paramref name="n"/>, rows 0 to <paramref name="n"/>-1 of <b>b</b> contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements <paramref name="n"/> to <paramref name="m"/>-1 in that column.</item>
        /// <item>If <paramref name="trans"/> == 'N' and <paramref name="m"/> &lt; <paramref name="n"/>, rows 0 to <paramref name="n"/>-1 of <b>b</b> contain the minimum norm solution vectors.</item>
        /// <item>If <paramref name="trans"/> == 'T' and <paramref name="m"/> &gt;= <paramref name="n"/>, rows 0 to <paramref name="m"/>-1 of <b>b</b> contain the minimum norm solution vectors.</item>
        /// <item>If <paramref name="trans"/> == 'T' and <paramref name="m"/> &lt;<paramref name="n"/>, rows 0 to <paramref name="m"/>-1 of <b>b</b> contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements <paramref name="m"/> to <paramref name="n"/>-1 in that column.</item>
        /// </list></param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <paramref name="ldb"/> &gt;= max(1,<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <param name="work">Work array. On exit, <b>work[0]</b> equals the optimal length the work array should be.
        /// The length of <b>work</b> should be at least max( 1, MN + max( MN, <paramref name="nrhs"/> ) ). For optimal performance,
        /// The length of <b>work</b> should be at least max( 1, MN + max( MN, <paramref name="nrhs"/> )*NB ), where MN = min(<paramref name="m"/>,<paramref name="n"/>) and NB is the optimum block size.
        /// </param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgels(char trans, int m, int n, int nrhs, double[] a, int lda, double[] b, int ldb, double[] work);

        /// <summary>
        /// <para>Computes the minimum norm solution to a real linear least squares problem:</para>
        /// <para>Minimize 2-norm(| B - A*X |)</para>
        /// <para>using the singular value decomposition (SVD) of A. A is an M by N matrix which may be rank-deficient.
        /// Several right hand side vectors B and solution vectors X can be handled in a single call. They are stored as the columns of the
        /// M by NRHS right hand side matrix B and the N by NRHS solution matrix X.
        /// The effective rank of A is determined by treating as zero those singular values which are less than rcond times the largest singular value.</para>
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns of the matrices <paramref name="b"/> and X. <paramref name="nrhs"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <paramref name="a"/>. On exit, the first min(<paramref name="m"/>,n) rows of <paramref name="a"/> 
        /// are overwritten with its right singular vectors, stored row wise.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <paramref name="lda"/> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b">Dimension <paramref name="ldb"/> by <paramref name="nrhs"/>. On entry, the <paramref name="m"/> by <paramref name="nrhs"/> right hand side matrix <paramref name="b"/>.
        /// On exit, <paramref name="b"/> is overwritten by the <paramref name="n"/> by <paramref name="nrhs"/> solution matrix X. If <paramref name="m"/> &gt;= <paramref name="n"/> and RANK = <paramref name="n"/>, the residual
        /// sum-of-squares for the solution in the i-th column is given	by the sum of squares of elements <paramref name="n"/>:<paramref name="m"/>-1 in that column.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <paramref name="ldb"/> &gt;= max(1,max(<paramref name="m"/>,<paramref name="n"/>)).</param>
        /// <param name="s">Dimension min(<paramref name="m"/>,<paramref name="n"/>) The singular values of <paramref name="a"/> in decreasing order. The condition number of <paramref name="a"/> 
        /// in the 2-norm = <paramref name="s"/>[0]/<paramref name="s"/>[min(<paramref name="m"/>,<paramref name="n"/>)-1].</param>
        /// <param name="rank">The effective rank of A.</param>
        /// <returns> <list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: the algorithm for computing the SVD failed to converge. i-th off-diagonal elements of an intermediate
        /// bidiagonal form did not converge to zero.</item>
        /// </list></returns>
        int Dgelss(int m, int n, int nrhs, double[] a, int lda, double[] b, int ldb, out double[] s, ref int rank);

        /// <summary>
        /// <para>Computes the minimum norm solution to a real linear least squares problem:</para>
        /// <para>Minimize 2-norm(| B - A*X |)</para>
        /// <para>using the singular value decomposition (SVD) of A. A is an M by N matrix which may be rank-deficient.
        /// Several right hand side vectors B and solution vectors X can be handled in a single call. They are stored as the columns of the
        /// M by NRHS right hand side matrix B and the N by NRHS solution matrix X.
        /// The effective rank of A is determined by treating as zero those singular values which are less than rcond times the largest singular value.</para>
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns of the matrices <paramref name="b"/> and X. <paramref name="nrhs"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <paramref name="a"/>. On exit, the first min(<paramref name="m"/>,n) rows of <paramref name="a"/> 
        /// are overwritten with its right singular vectors, stored row wise.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <paramref name="lda"/> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b">Dimension <paramref name="ldb"/> by <paramref name="nrhs"/>. On entry, the <paramref name="m"/> by <paramref name="nrhs"/> right hand side matrix <paramref name="b"/>.
        /// On exit, <paramref name="b"/> is overwritten by the <paramref name="n"/> by <paramref name="nrhs"/> solution matrix X. If <paramref name="m"/> &gt;= <paramref name="n"/> and RANK = <paramref name="n"/>, the residual
        /// sum-of-squares for the solution in the i-th column is given	by the sum of squares of elements <paramref name="n"/>:<paramref name="m"/>-1 in that column.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <paramref name="ldb"/> &gt;= max(1,max(<paramref name="m"/>,<paramref name="n"/>)).</param>
        /// <param name="s">Dimension min(<paramref name="m"/>,<paramref name="n"/>) The singular values of <paramref name="a"/> in decreasing order. The condition number of <paramref name="a"/> 
        /// in the 2-norm = S[0]/S[min(<paramref name="m"/>,<paramref name="n"/>)-1].</param>
        /// <param name="rcond"><b>rcond</b> is used to determine the effective rank of <paramref name="a"/>. Singular values <paramref name="s"/>[i] &lt;= <b>rcond</b>*<paramref name="s"/>[0] are treated as zero.
        /// If <b>rcond</b> &lt; 0, machine precision is used instead.</param>
        /// <param name="rank">The effective rank of A, i.e., the number of singular values which are greater than <paramref name="rcond"/>*<paramref name="s"/>[0].</param>
        /// <param name="work">Work array. On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use.
        /// The length of <b>work</b> should be at least 3*min(<paramref name="m"/>,<paramref name="n"/>) + max( 2*min(<paramref name="m"/>,<paramref name="n"/>), max(<paramref name="m"/>,<paramref name="n"/>), <paramref name="nrhs"/> ).
        /// </param>
        /// <returns> <list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: the algorithm for computing the SVD failed to converge. i-th off-diagonal elements of an intermediate
        /// bidiagonal form did not converge to zero.</item>
        /// </list></returns>
        int Dgelss(int m, int n, int nrhs, double[] a, int lda, double[] b, int ldb, out double[] s, double rcond, ref int rank, double[] work);

        /// <summary>
        /// <para>Computes the minimum-norm solution to a real linear least squares problem:</para>
        /// <para>minimize || A * X - B ||</para>
        /// <para>using a complete orthogonal factorization of A. A is an M by N
        /// matrix which may be rank-deficient. Several right hand side vectors b and solution vectors x can be
        /// handled in a single call; they are stored as the columns of the M by NRHS right hand side matrix B 
        /// and the N by NRHS solution matrix X. The routine first computes a QR factorization with column pivoting:
        /// <code>
        /// A * P = Q * [ R00 R01 ]
        ///             [ 0   R11 ]
        /// </code>
        /// with R11 defined as the largest leading submatrix whose estimated
        /// condition number is less than 1/RCOND. The order of R11, RANK,
        /// is the effective rank of A.</para>
        /// <para>Then, R22 is considered to be negligible, and R12 is annihilated
        /// by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
        /// <code>
        /// A * P = Q * [ T00 0 ] * Z
        ///             [ 0   0 ]
        /// </code>
        /// The minimum-norm solution is then
        /// <code>
        /// X = P * Z' [ inv(T00)*Q0'*B ]
        ///            [        0       ]
        /// </code>
        /// where Q0 consists of the first RANK columns of Q.</para>
        /// This routine is basically identical to the original xGELSX except three differences:
        /// <list>
        /// <item>The call to the subroutine xGEQPF has been substituted by the
        /// the call to the subroutine xGEQP3. This subroutine is a Blas-3
        /// version of the QR factorization with column pivoting.</item>
        /// <item>DenseMatrix B (the right hand side) is updated with Blas-3.</item>
        /// <item>The permutation of matrix B (the right hand side) is faster and more simple.</item>
        /// </list>
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of
        /// columns of matrices <paramref name="b"/> and X. <paramref name="nrhs"/> &gt;= 0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit, <b>a</b> has been overwritten by details of its complete orthogonal factorization.</param>
        /// <param name="lda"> The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b">dimension (<paramref name="ldb"/>,<paramref name="nrhs"/>)
        /// On entry, the <paramref name="m"/> by <paramref name="nrhs"/> right hand side matrix <paramref name="b"/>.
        /// On exit, the <paramref name="n"/> by <paramref name="nrhs"/> solution matrix X.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <b>ldb</b> &gt;= max(1,<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <param name="jpvt">On entry, if <b>jpvt[i] != 0</b>, the i-th column of <paramref name="a"/> is permuted
        /// to the front of <paramref name="a"/>P, otherwise column i is a free column.
        /// On exit, if <b>jpvt[i] = k</b>, then the i-th column of <paramref name="a"/>P
        /// was the k-th column of <paramref name="a"/>.</param>
        /// <param name="rcond"><b>rcond</b> is used to determine the effective rank of <paramref name="a"/>, which
        /// is defined as the order of the largest leading triangular
        /// submatrix R00 in the QR factorization with pivoting of <paramref name="a"/>,
        /// whose estimated condition number &lt; <b>1/rcond</b>.</param>
        /// <param name="rank"> The effective rank of <paramref name="a"/>, i.e., the order of the submatrix
        /// R00. This is the same as the order of the submatrix T00
        /// in the complete orthogonal factorization of <paramref name="a"/>.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgelsy(int m, int n, int nrhs, double[] a, int lda, double[] b, int ldb, int[] jpvt, double rcond, ref int rank);

        /// <summary>
        /// <para>Computes the minimum-norm solution to a real linear least squares problem:</para>
        /// <para>minimize || A * X - B ||</para>
        /// <para>using a complete orthogonal factorization of A. A is an M by N
        /// matrix which may be rank-deficient. Several right hand side vectors b and solution vectors x can be
        /// handled in a single call; they are stored as the columns of the M by NRHS right hand side matrix B 
        /// and the N by NRHS solution matrix X. The routine first computes a QR factorization with column pivoting:
        /// <code>
        /// A * P = Q * [ R00 R01 ]
        ///             [ 0   R11 ]
        /// </code>
        /// with R11 defined as the largest leading submatrix whose estimated
        /// condition number is less than 1/RCOND. The order of R11, RANK,
        /// is the effective rank of A.</para>
        /// <para>Then, R22 is considered to be negligible, and R12 is annihilated
        /// by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
        /// <code>
        /// A * P = Q * [ T00 0 ] * Z
        ///             [ 0   0 ]
        /// </code>
        /// The minimum-norm solution is then
        /// <code>
        /// X = P * Z' [ inv(T00)*Q0'*B ]
        ///            [        0       ]
        /// </code>
        /// where Q0 consists of the first RANK columns of Q.</para>
        /// This routine is basically identical to the original xGELSX except three differences:
        /// <list>
        /// <item>The call to the subroutine xGEQPF has been substituted by the
        /// the call to the subroutine xGEQP3. This subroutine is a Blas-3
        /// version of the QR factorization with column pivoting.</item>
        /// <item>DenseMatrix B (the right hand side) is updated with Blas-3.</item>
        /// <item>The permutation of matrix B (the right hand side) is faster and more simple.</item>
        /// </list>
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of
        /// columns of matrices <paramref name="b"/> and X. <paramref name="nrhs"/> &gt;= 0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit, <b>a</b> has been overwritten by details of its complete orthogonal factorization.</param>
        /// <param name="lda"> The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b">dimension (<paramref name="ldb"/>,<paramref name="nrhs"/>)
        /// On entry, the <paramref name="m"/> by <paramref name="nrhs"/> right hand side matrix <paramref name="b"/>.
        /// On exit, the <paramref name="n"/> by <paramref name="nrhs"/> solution matrix X.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <b>ldb</b> &gt;= max(1,<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <param name="jpvt">On entry, if <b>jpvt[i] != 0</b>, the i-th column of <paramref name="a"/> is permuted
        /// to the front of <paramref name="a"/>P, otherwise column i is a free column.
        /// On exit, if <b>jpvt[i] = k</b>, then the i-th column of <paramref name="a"/>P
        /// was the k-th column of <paramref name="a"/>.</param>
        /// <param name="rcond"><b>rcond</b> is used to determine the effective rank of <paramref name="a"/>, which
        /// is defined as the order of the largest leading triangular
        /// submatrix R00 in the QR factorization with pivoting of <paramref name="a"/>,
        /// whose estimated condition number &lt; <b>1/rcond</b>.</param>
        /// <param name="rank"> The effective rank of <paramref name="a"/>, i.e., the order of the submatrix
        /// R00. This is the same as the order of the submatrix T00
        /// in the complete orthogonal factorization of <paramref name="a"/>.</param>
        /// <param name="work">Work array. On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use.
        /// The length of <b>work</b> should be at least max( MN+2*<paramref name="n"/>+NB*(<paramref name="n"/>+1), 2*MN+NB*<paramref name="nrhs"/> ) 
        /// where MN = min( <paramref name="m"/>, <paramref name="n"/> ) and NB is an upper bound on the blocksize returned by ILAENV
        /// for the routines DGEQP3, DTZRZF, DTZRQF, DORMQR, and DORMRZ.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgelsy(int m, int n, int nrhs, double[] a, int lda, double[] b, int ldb, int[] jpvt, double rcond, ref int rank, double[] work);

        /// <summary>
        /// Computes a QR factorization with column pivoting of a matrix A: A*P = Q*R using Level 3 BLAS.
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit, the upper triangle of the array contains the
        /// min(<paramref name="m"/>,<paramref name="n"/>) by <paramref name="n"/> upper trapezoidal matrix R. The elements below
        /// the diagonal, together with the array <paramref name="tau"/>, represent the
        /// unitary matrix Q as a product of min(<paramref name="m"/>,<paramref name="n"/>) elementary
        /// reflectors.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="jpvt">Dimension <paramref name="n"/>.
        /// On entry, if <b>jpvt[j] != 0</b>, the j-th column of <paramref name="a"/> is permuted
        /// to the front of <paramref name="a"/>*P (a leading column); if <b>jpvt[j] == 0</b>,
        /// the j-th column of <paramref name="a"/> is a free column.
        /// On exit, if <b>jpvt[j] == k</b>, then the j-th column of <paramref name="a"/>*P was the
        /// the k-th column of <paramref name="a"/>.</param>
        /// <param name="tau">Dimension min(<paramref name="m"/>,<paramref name="n"/>). The scalar factors of the elementary reflectors.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgeqp3(int m, int n, double[] a, int lda, int[] jpvt, double[] tau);

        /// <summary>
        /// Computes a QR factorization with column pivoting of a matrix A: A*P = Q*R using Level 3 BLAS.
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit, the upper triangle of the array contains the
        /// min(<paramref name="m"/>,<paramref name="n"/>) by <paramref name="n"/> upper trapezoidal matrix R. The elements below
        /// the diagonal, together with the array <paramref name="tau"/>, represent the
        /// unitary matrix Q as a product of min(<paramref name="m"/>,<paramref name="n"/>) elementary
        /// reflectors.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="jpvt">Dimension <paramref name="n"/>.
        /// On entry, if <b>jpvt[j] != 0</b>, the j-th column of <paramref name="a"/> is permuted
        /// to the front of <paramref name="a"/>*P (a leading column); if <b>jpvt[j] == 0</b>,
        /// the j-th column of <paramref name="a"/> is a free column.
        /// On exit, if <b>jpvt[j] == k</b>, then the j-th column of <paramref name="a"/>*P was the
        /// the k-th column of <paramref name="a"/>.</param>
        /// <param name="tau">Dimension min(<paramref name="m"/>,<paramref name="n"/>). The scalar factors of the elementary reflectors.</param>
        /// <param name="work">Work array. On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use. 
        /// The length of <b>work</b> should be at least3*<paramref name="n"/>+1.
        /// For optimal performance the length of <b>work</b> should be at least2*<paramref name="n"/>+( <paramref name="n"/>+1 )*NB, where NB is the optimal blocksize.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgeqp3(int m, int n, double[] a, int lda, int[] jpvt, double[] tau, double[] work);

        /// <summary>
        /// Computes the singular value decomposition (SVD) of a complex
        /// M by N matrix A, optionally computing the left and/or right singular
        /// vectors. The SVD is written A = U * SIGMA * conjugate-transpose(V)
        /// where SIGMA is an M by N matrix which is zero except for its
        /// min(m,n) diagonal elements, U is an M by M unitary matrix, and	V is an N-by-N unitary matrix. The diagonal elements of SIGMA
        /// are the singular values of A. They are real and non-negative, and
        /// are returned in descending order. The first min(m,n) columns of
        /// U and V are the left and right singular vectors of A.
        /// Note that the routine returns V**T, not V.
        /// </summary>
        /// <param name="jobu">
        /// Specifies options for computing all or part of the matrix <paramref name="u"/>:
        /// <list>
        /// <item> = '<paramref name="a"/>': all <paramref name="m"/> columns of <paramref name="u"/> are returned in array <paramref name="u"/>.</item>
        /// <item> = '<paramref name="s"/>': the first min(m,n) columns of <paramref name="u"/> (the left singular vectors) are returned in the array <paramref name="u"/>.</item>
        /// <item> = 'O': the first min(m,n) columns of <paramref name="u"/> (the left singular vectors) are overwritten on the array <paramref name="a"/>.</item>
        /// <item> = '<paramref name="n"/>': no columns of <paramref name="u"/> (no left singular vectors) are computed. </item>
        /// </list>
        /// </param>
        /// <param name="jobvt">
        /// Specifies options for computing all or part of the matrix V**T:
        /// <list>
        /// <item>'<paramref name="a"/>': all <paramref name="n"/> rows of V**T are returned in the array <paramref name="vt"/>.</item>
        /// <item>'<paramref name="s"/>': the first min(m,n) rows of V**T (the right singular vectors) are returned in the array <paramref name="vt"/>.</item>
        /// <item>'O': the first min(m,n) rows of V**T (the right singular vectors) are overwritten on the array <paramref name="a"/>.</item>
        /// <item>'<paramref name="n"/>': no rows of V**T (no right singular vectors) are computed.</item>
        /// </list>
        /// </param>
        /// <param name="m">The number of rows of the input matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the input matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit:<list>
        /// <item>if <paramref name="jobu"/> == 'O', <b>a</b> is overwritten with the first min(m,n) columns of <paramref name="u"/> (the left singular vectors, stored column wise).</item>
        /// <item>if <paramref name="jobvt"/> == 'O', <b>a</b> is overwritten with the first min(m,n)rows of V**T (the right singular vectors, stored row wise).</item>
        /// <item>if <paramref name="jobu"/> != 'O' and <paramref name="jobvt"/> != 'O', the contents of <b>a</b> are destroyed.</item>
        /// </list></param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <param name="s">Dimension min(<paramref name="m"/>,<paramref name="n"/>). The singular values of <paramref name="a"/>, sorted so that <b>s[i] &gt;= s[i+1]</b>.</param>
        /// <param name="u">Dimension <paramref name="ldu"/> by UCOL, <paramref name="ldu"/> by <paramref name="m"/> if <paramref name="jobu"/> = '<paramref name="a"/>', or <paramref name="ldu"/> by min(<paramref name="m"/>,<paramref name="n"/>) if <paramref name="jobu"/> = '<paramref name="s"/>'.
        /// <list>
        /// <item>If <paramref name="jobu"/> == '<paramref name="a"/>', <b>u</b> contains the <paramref name="m"/> by <paramref name="m"/> orthogonal matrix <b>u</b>.</item>
        /// <item>If <paramref name="jobu"/> == '<paramref name="s"/>', <b>u</b> contains the first min(m,n) columns of <b>u</b> (the left singular vectors, stored column wise).</item>
        /// <item>if <paramref name="jobu"/> == '<paramref name="n"/>' or 'O', <b>u</b> is not referenced.</item>
        /// </list></param>
        /// <param name="ldu">The leading dimension of the array <paramref name="u"/>. <b>ldu &gt;= 1</b>. If <paramref name="jobu"/> = '<paramref name="s"/>' or '<paramref name="a"/>', <b>ldu &gt;= <paramref name="m"/></b>.</param>
        /// <param name="vt">Dimension <paramref name="ldvt"/> by <paramref name="n"/>.
        /// <list>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="a"/>', <b>vt</b> contains the <paramref name="n"/> by <paramref name="n"/> orthogonal matrix V**T.</item>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="s"/>', <b>vt</b> contains the first min(m,n) rows of V**T (the right singular vectors, stored row wise).</item>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="n"/>' or 'O', <b>vt</b> is not referenced.</item>
        /// </list></param>
        /// <param name="ldvt">The leading dimension of the array <paramref name="vt"/>. <b>ldvt &gt;= 1</b>. If <paramref name="jobvt"/> == '<paramref name="a"/>', 
        /// <b>ldvt &gt;= <paramref name="n"/></b>; if <paramref name="jobvt"/> == '<paramref name="s"/>', <b><paramref name="ldvt"/> &gt;= min(<paramref name="m"/>,<paramref name="n"/>)</b>.</param>
        /// <returns><list>
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: if DBDSQR did not converge, specifies how many
        /// superdiagonals of an intermediate bidiagonal form B
        /// did not converge to zero.</item>
        /// </list></returns>
        /// <remarks><paramref name="jobvt"/> and <paramref name="jobu"/> cannot both be 'O'.</remarks>
        int Dgesvd(char jobu, char jobvt, int m, int n, double[] a, int lda, double[] s, double[] u, int ldu, double[] vt, int ldvt);

        /// <summary>
        /// Computes the singular value decomposition (SVD) of a complex
        /// M by N matrix A, optionally computing the left and/or right singular
        /// vectors. The SVD is written A = U * SIGMA * conjugate-transpose(V)
        /// where SIGMA is an M by N matrix which is zero except for its
        /// min(m,n) diagonal elements, U is an M by M unitary matrix, and V is an N-by-N unitary matrix. The diagonal elements of SIGMA
        /// are the singular values of A. They are real and non-negative, and
        /// are returned in descending order. The first min(m,n) columns of
        /// U and V are the left and right singular vectors of A.
        /// Note that the routine returns V**T, not V.
        /// </summary>
        /// <param name="jobu">
        /// Specifies options for computing all or part of the matrix <paramref name="u"/>:
        /// <list>
        /// <item> = '<paramref name="a"/>': all <paramref name="m"/> columns of <paramref name="u"/> are returned in array <paramref name="u"/>.</item>
        /// <item> = '<paramref name="s"/>': the first min(m,n) columns of <paramref name="u"/> (the left singular vectors) are returned in the array <paramref name="u"/>.</item>
        /// <item> = 'O': the first min(m,n) columns of <paramref name="u"/> (the left singular vectors) are overwritten on the array <paramref name="a"/>.</item>
        /// <item> = '<paramref name="n"/>': no columns of <paramref name="u"/> (no left singular vectors) are computed. </item>
        /// </list>
        /// </param>
        /// <param name="jobvt">
        /// Specifies options for computing all or part of the matrix V**T:
        /// <list>
        /// <item>'<paramref name="a"/>': all <paramref name="n"/> rows of V**T are returned in the array <paramref name="vt"/>.</item>
        /// <item>'<paramref name="s"/>': the first min(m,n) rows of V**T (the right singular vectors) are returned in the array <paramref name="vt"/>.</item>
        /// <item>'O': the first min(m,n) rows of V**T (the right singular vectors) are overwritten on the array <paramref name="a"/>.</item>
        /// <item>'<paramref name="n"/>': no rows of V**T (no right singular vectors) are computed.</item>
        /// </list>
        /// </param>
        /// <param name="m">The number of rows of the input matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the input matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit:<list>
        /// <item>if <paramref name="jobu"/> == 'O', <b>a</b> is overwritten with the first min(m,n) columns of <paramref name="u"/> (the left singular vectors, stored column wise).</item>
        /// <item>if <paramref name="jobvt"/> == 'O', <b>a</b> is overwritten with the first min(m,n)rows of V**T (the right singular vectors, stored row wise).</item>
        /// <item>if <paramref name="jobu"/> != 'O' and <paramref name="jobvt"/> != 'O', the contents of <b>a</b> are destroyed.</item>
        /// </list></param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <param name="s">Dimension min(<paramref name="m"/>,<paramref name="n"/>). The singular values of <paramref name="a"/>, sorted so that <b>s[i] &gt;= s[i+1]</b>.</param>
        /// <param name="u">Dimension <paramref name="ldu"/> by UCOL, <paramref name="ldu"/> by <paramref name="m"/> if <paramref name="jobu"/> = '<paramref name="a"/>', or <paramref name="ldu"/> by min(<paramref name="m"/>,<paramref name="n"/>) if <paramref name="jobu"/> = '<paramref name="s"/>'.
        /// <list>
        /// <item>If <paramref name="jobu"/> == '<paramref name="a"/>', <b>u</b> contains the <paramref name="m"/> by <paramref name="m"/> orthogonal matrix <b>u</b>.</item>
        /// <item>If <paramref name="jobu"/> == '<paramref name="s"/>', <b>u</b> contains the first min(m,n) columns of <b>u</b> (the left singular vectors, stored column wise).</item>
        /// <item>if <paramref name="jobu"/> == '<paramref name="n"/>' or 'O', <b>u</b> is not referenced.</item>
        /// </list></param>
        /// <param name="ldu">The leading dimension of the array <paramref name="u"/>. <b>ldu &gt;= 1</b>. If <paramref name="jobu"/> = '<paramref name="s"/>' or '<paramref name="a"/>', <b>ldu &gt;= <paramref name="m"/></b>.</param>
        /// <param name="vt">Dimension <paramref name="ldvt"/> by <paramref name="n"/>.
        /// <list>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="a"/>', <b>vt</b> contains the <paramref name="n"/> by <paramref name="n"/> orthogonal matrix V**T.</item>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="s"/>', <b>vt</b> contains the first min(m,n) rows of V**T (the right singular vectors, stored row wise).</item>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="n"/>' or 'O', <b>vt</b> is not referenced.</item>
        /// </list></param>
        /// <param name="ldvt">The leading dimension of the array <paramref name="vt"/>. <b>ldvt &gt;= 1</b>. If <paramref name="jobvt"/> == '<paramref name="a"/>', 
        /// <b>ldvt &gt;= <paramref name="n"/></b>; if <paramref name="jobvt"/> == '<paramref name="s"/>', <b><paramref name="ldvt"/> &gt;= min(<paramref name="m"/>,<paramref name="n"/>)</b>.</param>
        /// <param name="work">Work array. The length of <b>work</b> should be at least max(3*min(<paramref name="m"/>,<paramref name="n"/>)+max(<paramref name="m"/>,<paramref name="n"/>),5*min(<paramref name="m"/>,<paramref name="n"/>)).
        /// On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use. 
        /// If the method returns a value greater than 0, <b>work[1:min(<paramref name="m"/>,<paramref name="n"/>)-1]</b> contains the unconverged superdiagonal elements of an upper bidiagonal
        /// matrix B whose diagonal is in <paramref name="s"/> (not necessarily sorted).
        /// B satisfies <paramref name="a"/> = <paramref name="u"/> * B * <paramref name="vt"/>, so it has the same singular
        /// values as <paramref name="a"/>, and singular vectors related by <paramref name="u"/> and <paramref name="vt"/>.</param>
        /// <returns><list>
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: if DBDSQR did not converge, specifies how many
        /// superdiagonals of an intermediate bidiagonal form B
        /// did not converge to zero.See the description of <paramref name="work"/>
        /// above for details.</item>
        /// </list></returns>
        /// <remarks><paramref name="jobvt"/> and <paramref name="jobu"/> cannot both be 'O'.</remarks>
        int Dgesvd(char jobu, char jobvt, int m, int n, double[] a, int lda, double[] s, double[] u, int ldu, double[] vt, int ldvt, double[] work);

        /// <summary>
        /// Computes an LU factorization of a general M by N matrix A
        /// using partial pivoting with row interchanges. 
        /// The factorization has the form A = P * L * U
        /// where P is a permutation matrix, L is lower triangular with unit
        /// diagonal elements (lower trapezoidal if m &gt; n), and U is upper
        /// triangular (upper trapezoidal if m &lt; n).
        /// This is the right-looking Level 3 BLAS version of the algorithm.
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimensions <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/>-by-<paramref name="n"/> matrix to be factored.
        /// On exit, the factors L and U from the factorization
        /// A = P*L*U. The unit diagonal elements of L are not stored.</param>
        /// <param name="lda"> The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <param name="ipiv">Dimension min(<paramref name="m"/>,<paramref name="n"/>).
        /// The pivot indices. for 0 &lt;= i &lt;= min(<paramref name="m"/>,<paramref name="n"/>), row i of the matrix was interchanged with row <b>ipiv[i]</b>.</param>
        /// <returns>
        /// <list> 
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: U[i,i] is exactly zero. The factorization has been completed, but the factor U is exactly
        ///	singular, and division by zero will occur if it is used
        /// to solve a system of equations.</item>
        /// </list></returns>
        int Dgetrf(int m, int n, double[] a, int lda, out int[] ipiv);

        /// <summary>
        /// Computes the inverse of a matrix using the LU factorization
        /// computed by DGETRF. This method inverts U and then computes 
        /// inv(A) by solving the system inv(A)*L = inv(U) for inv(A).
        /// </summary>
        /// <param name="n"> The order of the matrix <parmaref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <parmaref name="lda"/> by <parmaref name="n"/>. On entry, the factors L and U from the factorization
        /// A = P*L*U as computed by DGETRF.
        /// On exit, the inverse of the original matrix <b>a</b>.</param>
        /// <param name="lda">The leading dimension of the array <parmaref name="a"/>. <b>lda</b> &gt;= max(1,<parmaref name="n"/>).</param>
        /// <param name="ipiv">Dimension <parmaref name="n"/>. The pivot indices from DGETRF; for 0 &lt;= i &lt;= <parmaref name="n"/>-1, row i of the
        /// matrix was interchanged with row <b>ipiv[i]</b>.</param>
        /// <returns><list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value</item>
        /// <item>&gt; 0: U[i,i] is exactly zero; the matrix is
        /// singular and its inverse could not be computed.</item>
        /// </list></returns>
        int Dgetri(int n, double[] a, int lda, int[] ipiv);

        /// <summary>
        /// Computes the inverse of a matrix using the LU factorization
        /// computed by DGETRF. This method inverts U and then computes 
        /// inv(A) by solving the system inv(A)*L = inv(U) for inv(A).
        /// </summary>
        /// <param name="n"> The order of the matrix <parmaref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <parmaref name="lda"/> by <parmaref name="n"/>. On entry, the factors L and U from the factorization
        /// A = P*L*U as computed by DGETRF.
        /// On exit, the inverse of the original matrix <b>a</b>.</param>
        /// <param name="lda">The leading dimension of the array <parmaref name="a"/>. <b>lda</b> &gt;= max(1,<parmaref name="n"/>).</param>
        /// <param name="ipiv">Dimension <parmaref name="n"/>. The pivot indices from DGETRF; for 0 &lt;= i &lt;= <parmaref name="n"/>-1, row i of the
        /// matrix was interchanged with row <b>ipiv[i]</b>.</param>
        /// <param name="work">Work array. The length of <b>work</b> should be at least <paramref name="n"/>*NB,
        /// where NB is the optimal blocksize returned by ILAENV.
        /// On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use.</param>
        /// <returns><list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value</item>
        /// <item>&gt; 0: U[i,i] is exactly zero; the matrix is
        /// singular and its inverse could not be computed.</item>
        /// </list></returns>
        int Dgetri(int n, double[] a, int lda, int[] ipiv, double[] work);

        /// <summary>
        /// Solves a system of linear equations A * X = B or A' * X = B
        /// with a general N by N matrix A using the LU factorization computed by DGETRF.
        /// </summary>
        /// <param name="trans"> Specifies the form of the system of equations:
        /// <list>
        /// <item>'N': A * X = B (No transpose)</item>
        /// <item>'T': A'* X = B (Transpose)</item>
        /// <item>'C': A'* X = B (Conjugate transpose = Transpose)</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns
        /// of the matrix <paramref name="b"/>. <b>nrhs &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// The factors L and U from the factorization A = P*L*U as computed by DGETRF.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <param name="ipiv">Dimension <parmaref name="n"/>. The pivot indices from DGETRF; for 0 &lt;= i &lt;= <parmaref name="n"/>-1, row i of the
        /// matrix was interchanged with row <b>ipiv[i]</b>.</param>
        /// <param name="b">Dimension <parmaref name="ldb"/> by <parmaref name="nrhs"/>.
        /// On entry, the right hand side matrix <b>b</b>.
        /// On exit, the solution matrix X.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <b>ldb</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgetrs(char trans, int n, int nrhs, double[] a, int lda, int[] ipiv, double[] b, int ldb);

        /// <summary>
        /// Generates an M by N real matrix Q with orthonormal columns,
        /// which is defined as the first N columns of a product of K elementary
        /// reflectors of order M: <code>Q = H(0) H(1) . . . H(k-1)</code>
        /// as returned by DGEQRF.
        /// </summary>
        /// <param name="m">The number of rows of the matrix Q. <b>m &gt;= 0</b>.</param>
        /// <param name="n">The number of columns of the matrix Q. <b>m &gt;= n &gt;= 0</b></param>
        /// <param name="k">The number of elementary reflectors whose product defines the matrix Q. <b>n &gt;= k &gt;= 0</b>.</param>
        /// <param name="a">Dimension LDA by N.
        /// On entry, the i-th column must contain the vector which defines the elementary reflector H(0), for i = 0,2,...,k-1, as
        /// returned by DGEQRF in the first k columns of its array argument <b>a</b>.
        /// On exit, the <parmaref name="m"/> by <parmaref name="n"/> matrix Q.</param>
        /// <param name="lda">The first dimension of the array <parmaref name="a"/> . <b>lda</b> &gt;= max(1,<parmaref name="m"/>).</param>
        /// <param name="tau">Dimension K. <b>tau[i]</b> must contain the scalar factor of the elementary
        /// reflector H(i), as returned by DGEQRF.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dorgqr(int m, int n, int k, double[] a, int lda, double[] tau);

        /// <summary>
        /// Generates an M by N real matrix Q with orthonormal columns,
        /// which is defined as the first N columns of a product of K elementary
        /// reflectors of order M: <code>Q = H(0) H(1) . . . H(k-1)</code>
        /// as returned by DGEQRF.
        /// </summary>
        /// <param name="m">The number of rows of the matrix Q. <b>m &gt;= 0</b>.</param>
        /// <param name="n">The number of columns of the matrix Q. <b>m &gt;= n &gt;= 0</b></param>
        /// <param name="k">The number of elementary reflectors whose product defines the matrix Q. <b>n &gt;= k &gt;= 0</b>.</param>
        /// <param name="a">Dimension LDA by N.
        /// On entry, the i-th column must contain the vector which defines the elementary reflector H(0), for i = 0,2,...,k-1, as
        /// returned by DGEQRF in the first k columns of its array argument <b>a</b>.
        /// On exit, the <parmaref name="m"/> by <parmaref name="n"/> matrix Q.</param>
        /// <param name="lda">The first dimension of the array <parmaref name="a"/> . <b>lda</b> &gt;= max(1,<parmaref name="m"/>).</param>
        /// <param name="tau">Dimension K. <b>tau[i]</b> must contain the scalar factor of the elementary
        /// reflector H(i), as returned by DGEQRF.</param>
        /// <param name="work">Work array. The length of <b>work</b> should be at least <paramref name="n"/>*NB,
        /// where NB is the optimal blocksize returned by ILAENV.
        /// On exit, <b>work[0]</b> returns the optimal length of <b>work</b> to use.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dorgqr(int m, int n, int k, double[] a, int lda, double[] tau, double[] work);

        /// <summary>
        /// <para>Overwrites the general real M-by-N matrix C with</para>
        /// <pre>             SIDE = 'L'    SIDE = 'R'</pre>
        /// <pre>TRANS = 'N':   Q * C         C * Q</pre>
        /// <pre>TRANS = 'T':   Q**T * C      C * Q**T</pre>																		
        /// <para>where Q is a real orthogonal matrix defined as the product of k elementary reflectors
        /// Q = H(0) H(1) . . . H(k-1) as returned by DGEQRF. Q is of order M if 
        /// SIDE == 'L' and of order N if SIDE == 'R'.</para>
        /// </summary>
        /// <param name="side"><list>
        /// <item>'L': apply Q or Q**T from the Left.</item>
        /// <item>'R': apply Q or Q**T from the Right.</item>
        /// </list></param>
        /// <param name="trans"><list>
        /// <item>'N': No transpose, apply Q;</item>
        /// <item>'T': Transpose, apply Q**T.</item>
        /// </list></param>
        /// <param name="m">The number of rows of the matrix <paramref name="c"/>. <b>m &gt;= 0</b>.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="c"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="k"> The number of elementary reflectors whose product defines the matrix Q.
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <paramref name="m"/> &gt;= <paramref name="k"/> &gt;= 0.</item>
        /// <item>If <paramref name="side"/> == 'R', <paramref name="n"/> &gt;= <paramref name="k"/> &gt;= 0.</item>
        /// </list></param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="k"/>.
        /// The i-th column must contain the vector which defines the
        /// elementary reflector H(i), for i = 0,2,...,<paramref name="k"/>-1, as returned by
        ///	DGEQRF in the first <paramref name="k"/> columns of its array argument <paramref name="a"/>.
        ///	<paramref name="a"/> is modified by the routine but restored on exit.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/> .
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <b>lda</b> &gt;= max(1,<paramref name="m"/>).</item>
        /// <item>If <paramref name="side"/> == 'R', <b>lda</b> &gt;= max(1,<paramref name="n"/>).</item>
        /// </list></param>
        /// <param name="tau">Dimension <paramref name="k"/>. <b>tau[i]</b> must contain the scalar factor of the elementary
        /// reflector H(i), as returned by DGEQRF.</param>
        /// <param name="c">dimension <paramref name="ldc"/> by <paramref name="n"/>. 
        /// On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>c</b>.
        /// On exit, <b>c</b> is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.</param>
        /// <param name="ldc">The leading dimension of the array <paramref name="c"/>.<b>ldc</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dormqr(char side, char trans, int m, int n, int k, double[] a, int lda, double[] tau, double[] c, int ldc);

        /// <summary>
        /// <para>Overwrites the general real M-by-N matrix C with</para>
        /// <pre>             SIDE = 'L'    SIDE = 'R'</pre>
        /// <pre>TRANS = 'N':   Q * C         C * Q</pre>
        /// <pre>TRANS = 'T':   Q**T * C      C * Q**T</pre>																		
        /// <para>where Q is a real orthogonal matrix defined as the product of k elementary reflectors
        /// Q = H(0) H(1) . . . H(k-1) as returned by DGEQRF. Q is of order M if 
        /// SIDE == 'L' and of order N if SIDE == 'R'.</para>
        /// </summary>
        /// <param name="side"><list>
        /// <item>'L': apply Q or Q**T from the Left.</item>
        /// <item>'R': apply Q or Q**T from the Right.</item>
        /// </list></param>
        /// <param name="trans"><list>
        /// <item>'N': No transpose, apply Q;</item>
        /// <item>'T': Transpose, apply Q**T.</item>
        /// </list></param>
        /// <param name="m">The number of rows of the matrix <paramref name="c"/>. <b>m &gt;= 0</b>.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="c"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="k"> The number of elementary reflectors whose product defines the matrix Q.
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <paramref name="m"/> &gt;= <paramref name="k"/> &gt;= 0.</item>
        /// <item>If <paramref name="side"/> == 'R', <paramref name="n"/> &gt;= <paramref name="k"/> &gt;= 0.</item>
        /// </list></param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="k"/>.
        /// The i-th column must contain the vector which defines the
        /// elementary reflector H(i), for i = 0,2,...,<paramref name="k"/>-1, as returned by
        ///	DGEQRF in the first <paramref name="k"/> columns of its array argument <paramref name="a"/>.
        ///	<paramref name="a"/> is modified by the routine but restored on exit.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/> .
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <b>lda</b> &gt;= max(1,<paramref name="m"/>).</item>
        /// <item>If <paramref name="side"/> == 'R', <b>lda</b> &gt;= max(1,<paramref name="n"/>).</item>
        /// </list></param>
        /// <param name="tau">Dimension <paramref name="k"/>. <b>tau[i]</b> must contain the scalar factor of the elementary
        /// reflector H(i), as returned by DGEQRF.</param>
        /// <param name="c">dimension <paramref name="ldc"/> by <paramref name="n"/>. 
        /// On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>c</b>.
        /// On exit, <b>c</b> is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.</param>
        /// <param name="ldc">The leading dimension of the array <paramref name="c"/>.<b>ldc</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="work">Work array. The length of <b>work</b> should be:
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <paramref name="n"/> * NB.</item>
        /// <item>if <paramref name="side"/> == 'R', <paramref name="m"/> * NB.</item>
        /// </list>
        /// where NB is the optimal blocksize returned by ILAENV.
        /// On exit, <b>work[0]</b> returns the optimal length of <b>work</b> to use.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dormqr(char side, char trans, int m, int n, int k, double[] a, int lda, double[] tau, double[] c, int ldc, double[] work);

        /// <summary>
        /// Computes the Cholesky factorization of a real symmetric
        /// positive definite matrix A.
        /// The factorization has the form 
        /// <code>
        /// A = U**T * U, if UPLO == 'U' 
        /// A = L * L**T, if UPLO == 'L'
        /// </code>
        /// where U is an upper triangular matrix and L is lower triangular.
        /// This is the block version of the algorithm, calling Level 3 BLAS.
        /// </summary>
        /// <param name="uplo"><list>
        /// <item>'U': Upper triangle of A is stored.</item>
        /// <item>'L': Lower triangle of A is stored.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// On entry, the symmetric matrix <b>a</b>. 
        /// If <paramref name="uplo"/> == 'U', the leading <paramref name="n"/> by <paramref name="n"/> upper triangular part of <b>a</b>
        /// contains the upper triangular part of the matrix <b>a</b>, and the strictly lower triangular part of <b>a</b> is not referenced.
        /// If <paramref name="uplo"/> == 'L', the leading <paramref name="n"/> by <paramref name="n"/> lower triangular part of <b>a</b>
        /// contains the lower triangular part of the matrix <b>a</b>, and the strictly upper triangular part of <b>a</b> is not referenced.
        /// On exit, the factor U or L from the Cholesky factorization <b>A = U**T*U or A = L*L**T</b>.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <returns><list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value</item>
        /// <item>&gt; 0: the leading minor of order i is not positive definite, and the factorization could not be completed.</item>
        /// </list></returns>
        int Dpotrf(char uplo, int n, double[] a, int lda);

        /// <summary>
        /// Computes the inverse of a real symmetric positive definite
        /// matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
        /// computed by DPOTRF.
        /// </summary>
        /// <param name="uplo"><list>
        /// <item>'U': Upper triangle of A is stored.</item>
        /// <item>'L': Lower triangle of A is stored.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// On entry, the triangular factor U or L from the Cholesky
        /// factorization A = U**T*U or A = L*L**T, as computed by DPOTRF.
        /// On exit, the upper or lower triangle of the (symmetric) inverse of <b>a</b>, 
        /// overwriting the input factor U or L.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <returns><list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: the [i,i] element of the factor U or L is zero, and the inverse could not be computed.</item>
        /// </list></returns>
        int Dpotri(char uplo, int n, double[] a, int lda);

        /// <summary>
        /// Solves a system of linear equations A*X = B with a symmetric
        /// positive definite matrix A using the Cholesky factorization
        /// <code>A = U**T*U or A = L*L**T</code> computed by DPOTRF.
        /// </summary>
        /// <param name="uplo"><list>
        /// <item>'U': Upper triangle of A is stored.</item>
        /// <item>'L': Lower triangle of A is stored.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns
        /// of the matrix <paramref name="b"/>. <b>nrhs &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// The triangular factor U or L from the Cholesky factorization
        /// A = U**T*U or A = L*L**T, as computed by DPOTRF.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <param name="b">Dimension <paramref name="ldb"/> by <paramref name="nrhs"/>.
        /// On entry, the right hand side matrix <b>b</b>.
        /// On exit, the solution matrix X.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <b>ldb</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dpotrs(char uplo, int n, int nrhs, double[] a, int lda, double[] b, int ldb);

        /// <summary>
        /// Computes the least squares solution to an over-determined system of linear equations, A*X=B, or the minimum norm
        /// solution of an under-determined system, where A is a general rectangular matrix of full rank, using a QR or LQ factorization
        /// of A or under determined linear system with full rank matrix.
        /// </summary>
        /// <param name="trans">If <b>trans</b> == 'N': the linear system involves <paramref name="a"/>. If <b>trans</b> == 'T': the linear system involves <paramref name="a"/>'.</param>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <b>m</b> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <b>n</b> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns of the matrix <paramref name="b"/>. <b>nrhs</b> &gt;=0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix. On exit, if <paramref name="m"/> >= <paramref name="n"/>, 
        /// <paramref name="a"/> is overwritten by details of its QR factorization as returned by Sgeqrf. Otherwise, <paramref name="a"/> is 
        /// overwritten by details of its LQ factorization as returned by Sgelqf.</param>
        /// <param name="lda">The leading dimension of <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b"> Dimension (<paramref name="ldb"/>,<paramref name="nrhs"/>). On entry, the matrix <b>b</b> of right hand side vectors, stored column wise. 
        /// b is <paramref name="m"/> by <paramref name="nrhs"/> if <paramref name="trans"/> == 'N', or <paramref name="n"/> by <paramref name="n"/> by <paramref name="nrhs"/>
        /// if <paramref name="trans"/> == 'T'. On exit, <b>b</b> is overwritten by the solution vectors, stored column wise: 
        /// <list>
        /// <item>If <paramref name="trans"/> == 'N' and <paramref name="m"/> &gt;= <paramref name="n"/>, rows 0 to <paramref name="n"/>-1 of <b>b</b> contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements <paramref name="n"/> to <paramref name="m"/>-1 in that column.</item>
        /// <item>If <paramref name="trans"/> == 'N' and <paramref name="m"/> &lt; <paramref name="n"/>, rows 0 to <paramref name="n"/>-1 of <b>b</b> contain the minimum norm solution vectors.</item>
        /// <item>If <paramref name="trans"/> == 'T' and <paramref name="m"/> &gt;= <paramref name="n"/>, rows 0 to <paramref name="m"/>-1 of <b>b</b> contain the minimum norm solution vectors.</item>
        /// <item>If <paramref name="trans"/> == 'T' and <paramref name="m"/> &lt;<paramref name="n"/>, rows 0 to <paramref name="m"/>-1 of <b>b</b> contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements <paramref name="m"/> to <paramref name="n"/>-1 in that column.</item>
        /// </list></param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <paramref name="ldb"/> &gt;= max(1,<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgels(char trans, int m, int n, int nrhs, Complex[] a, int lda, Complex[] b, int ldb);

        /// <summary>
        /// Computes the least squares solution to an over-determined system of linear equations, A*X=B, or the minimum norm
        /// solution of an under-determined system, where A is a general rectangular matrix of full rank, using a QR or LQ factorization
        /// of A or under determined linear system with full rank matrix.
        /// </summary>
        /// <param name="trans">If <b>trans</b> == 'N': the linear system involves <paramref name="a"/>. If <b>trans</b> == 'T': the linear system involves <paramref name="a"/>'.</param>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <b>m</b> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <b>n</b> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns of the matrix <paramref name="b"/>. <b>nrhs</b> &gt;=0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix. On exit, if <paramref name="m"/> >= <paramref name="n"/>, 
        /// <paramref name="a"/> is overwritten by details of its QR factorization as returned by Sgeqrf. Otherwise, <paramref name="a"/> is 
        /// overwritten by details of its LQ factorization as returned by Sgelqf.</param>
        /// <param name="lda">The leading dimension of <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b"> Dimension (<paramref name="ldb"/>,<paramref name="nrhs"/>). On entry, the matrix <b>b</b> of right hand side vectors, stored column wise. 
        /// b is <paramref name="m"/> by <paramref name="nrhs"/> if <paramref name="trans"/> == 'N', or <paramref name="n"/> by <paramref name="n"/> by <paramref name="nrhs"/>
        /// if <paramref name="trans"/> == 'T'. On exit, <b>b</b> is overwritten by the solution vectors, stored column wise: 
        /// <list>
        /// <item>If <paramref name="trans"/> == 'N' and <paramref name="m"/> &gt;= <paramref name="n"/>, rows 0 to <paramref name="n"/>-1 of <b>b</b> contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements <paramref name="n"/> to <paramref name="m"/>-1 in that column.</item>
        /// <item>If <paramref name="trans"/> == 'N' and <paramref name="m"/> &lt; <paramref name="n"/>, rows 0 to <paramref name="n"/>-1 of <b>b</b> contain the minimum norm solution vectors.</item>
        /// <item>If <paramref name="trans"/> == 'T' and <paramref name="m"/> &gt;= <paramref name="n"/>, rows 0 to <paramref name="m"/>-1 of <b>b</b> contain the minimum norm solution vectors.</item>
        /// <item>If <paramref name="trans"/> == 'T' and <paramref name="m"/> &lt;<paramref name="n"/>, rows 0 to <paramref name="m"/>-1 of <b>b</b> contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements <paramref name="m"/> to <paramref name="n"/>-1 in that column.</item>
        /// </list></param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <paramref name="ldb"/> &gt;= max(1,<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <param name="work">Work array. On exit, <b>work[0]</b> equals the optimal length the work array should be.
        /// The length of <b>work</b> should be at least max( 1, MN + max( MN, <paramref name="nrhs"/> ) ). For optimal performance,
        /// The length of <b>work</b> should be at least max( 1, MN + max( MN, <paramref name="nrhs"/> )*NB ), where MN = min(<paramref name="m"/>,<paramref name="n"/>) and NB is the optimum block size.
        /// </param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgels(char trans, int m, int n, int nrhs, Complex[] a, int lda, Complex[] b, int ldb, Complex[] work);

        /// <summary>
        /// <para>Computes the minimum norm solution to a real linear least squares problem:</para>
        /// <para>Minimize 2-norm(| B - A*X |)</para>
        /// <para>using the singular value decomposition (SVD) of A. A is an M by N matrix which may be rank-deficient.
        /// Several right hand side vectors B and solution vectors X can be handled in a single call. They are stored as the columns of the
        /// M by NRHS right hand side matrix B and the N by NRHS solution matrix X.
        /// The effective rank of A is determined by treating as zero those singular values which are less than rcond times the largest singular value.</para>
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns of the matrices <paramref name="b"/> and X. <paramref name="nrhs"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <paramref name="a"/>. On exit, the first min(<paramref name="m"/>,n) rows of <paramref name="a"/> 
        /// are overwritten with its right singular vectors, stored row wise.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <paramref name="lda"/> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b">Dimension <paramref name="ldb"/> by <paramref name="nrhs"/>. On entry, the <paramref name="m"/> by <paramref name="nrhs"/> right hand side matrix <paramref name="b"/>.
        /// On exit, <paramref name="b"/> is overwritten by the <paramref name="n"/> by <paramref name="nrhs"/> solution matrix X. If <paramref name="m"/> &gt;= <paramref name="n"/> and RANK = <paramref name="n"/>, the residual
        /// sum-of-squares for the solution in the i-th column is given	by the sum of squares of elements <paramref name="n"/>:<paramref name="m"/>-1 in that column.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <paramref name="ldb"/> &gt;= max(1,max(<paramref name="m"/>,<paramref name="n"/>)).</param>
        /// <param name="s">Dimension min(<paramref name="m"/>,<paramref name="n"/>) The singular values of <paramref name="a"/> in decreasing order. The condition number of <paramref name="a"/> 
        /// in the 2-norm = <paramref name="s"/>[0]/<paramref name="s"/>[min(<paramref name="m"/>,<paramref name="n"/>)-1].</param>
        /// <param name="rank">The effective rank of A.</param>
        /// <returns> <list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: the algorithm for computing the SVD failed to converge. i-th off-diagonal elements of an intermediate
        /// bidiagonal form did not converge to zero.</item>
        /// </list></returns>
        int Zgelss(int m, int n, int nrhs, Complex[] a, int lda, Complex[] b, int ldb, out double[] s, ref int rank);

        /// <summary>
        /// <para>Computes the minimum norm solution to a real linear least squares problem:</para>
        /// <para>Minimize 2-norm(| B - A*X |)</para>
        /// <para>using the singular value decomposition (SVD) of A. A is an M by N matrix which may be rank-deficient.
        /// Several right hand side vectors B and solution vectors X can be handled in a single call. They are stored as the columns of the
        /// M by NRHS right hand side matrix B and the N by NRHS solution matrix X.
        /// The effective rank of A is determined by treating as zero those singular values which are less than rcond times the largest singular value.</para>
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns of the matrices <paramref name="b"/> and X. <paramref name="nrhs"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <paramref name="a"/>. On exit, the first min(<paramref name="m"/>,n) rows of <paramref name="a"/> 
        /// are overwritten with its right singular vectors, stored row wise.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <paramref name="lda"/> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b">Dimension <paramref name="ldb"/> by <paramref name="nrhs"/>. On entry, the <paramref name="m"/> by <paramref name="nrhs"/> right hand side matrix <paramref name="b"/>.
        /// On exit, <paramref name="b"/> is overwritten by the <paramref name="n"/> by <paramref name="nrhs"/> solution matrix X. If <paramref name="m"/> &gt;= <paramref name="n"/> and RANK = <paramref name="n"/>, the residual
        /// sum-of-squares for the solution in the i-th column is given	by the sum of squares of elements <paramref name="n"/>:<paramref name="m"/>-1 in that column.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <paramref name="ldb"/> &gt;= max(1,max(<paramref name="m"/>,<paramref name="n"/>)).</param>
        /// <param name="s">Dimension min(<paramref name="m"/>,<paramref name="n"/>) The singular values of <paramref name="a"/> in decreasing order. The condition number of <paramref name="a"/> 
        /// in the 2-norm = S[0]/S[min(<paramref name="m"/>,<paramref name="n"/>)-1].</param>
        /// <param name="rcond"><b>rcond</b> is used to determine the effective rank of <paramref name="a"/>. Singular values <paramref name="s"/>[i] &lt;= <b>rcond</b>*<paramref name="s"/>[0] are treated as zero.
        /// If <b>rcond</b> &lt; 0, machine precision is used instead.</param>
        /// <param name="rank">The effective rank of A, i.e., the number of singular values which are greater than <paramref name="rcond"/>*<paramref name="s"/>[0].</param>
        /// <param name="work">Work array. On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use.
        /// The length of <b>work</b> should be at least 2*min(<paramref name="m"/>,<paramref name="n"/>) + max(<paramref name="m"/>,<paramref name="n"/>,<paramref name="nrhs"/>).
        /// </param>
        /// <param name="rwork">Dimension 5*min(<paramref name="m"/>,<paramref name="n"/>)</param>
        /// <returns> <list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: the algorithm for computing the SVD failed to converge. i-th off-diagonal elements of an intermediate
        /// bidiagonal form did not converge to zero.</item>
        /// </list></returns>
        int Zgelss(int m, int n, int nrhs, Complex[] a, int lda, Complex[] b, int ldb, out double[] s, double rcond, ref int rank, Complex[] work, double[] rwork);

        /// <summary>
        /// <para>Computes the minimum-norm solution to a real linear least squares problem:</para>
        /// <para>minimize || A * X - B ||</para>
        /// <para>using a complete orthogonal factorization of A. A is an M by N
        /// matrix which may be rank-deficient. Several right hand side vectors b and solution vectors x can be
        /// handled in a single call; they are stored as the columns of the M by NRHS right hand side matrix B 
        /// and the N by NRHS solution matrix X. The routine first computes a QR factorization with column pivoting:
        /// <code>
        /// A * P = Q * [ R00 R01 ]
        ///             [ 0   R11 ]
        ///</code>
        /// with R11 defined as the largest leading submatrix whose estimated
        /// condition number is less than 1/RCOND. The order of R11, RANK,
        /// is the effective rank of A.</para>
        /// <para>Then, R22 is considered to be negligible, and R12 is annihilated
        /// by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
        /// <code>
        /// A * P = Q * [ T00 0 ] * Z
        ///             [ 0   0 ]
        /// </code>
        /// The minimum-norm solution is then
        /// <code>
        /// X = P * Z' [ inv(T00)*Q0'*B ]
        ///            [        0       ]
        /// </code>
        /// where Q0 consists of the first RANK columns of Q.</para>
        /// This routine is basically identical to the original xGELSX except three differences:
        /// <list>
        /// <item>The call to the subroutine xGEQPF has been substituted by the
        /// the call to the subroutine xGEQP3. This subroutine is a Blas-3
        /// version of the QR factorization with column pivoting.</item>
        /// <item>DenseMatrix B (the right hand side) is updated with Blas-3.</item>
        /// <item>The permutation of matrix B (the right hand side) is faster and more simple.</item>
        /// </list>
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of
        /// columns of matrices <paramref name="b"/> and X. <paramref name="nrhs"/> &gt;= 0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit, <b>a</b> has been overwritten by details of its complete orthogonal factorization.</param>
        /// <param name="lda"> The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b">dimension (<paramref name="ldb"/>,<paramref name="nrhs"/>)
        /// On entry, the <paramref name="m"/> by <paramref name="nrhs"/> right hand side matrix <paramref name="b"/>.
        /// On exit, the <paramref name="n"/> by <paramref name="nrhs"/> solution matrix X.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <b>ldb</b> &gt;= max(1,<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <param name="jpvt">On entry, if <b>jpvt[i] != 0</b>, the i-th column of <paramref name="a"/> is permuted
        /// to the front of <paramref name="a"/>P, otherwise column i is a free column.
        /// On exit, if <b>jpvt[i] = k</b>, then the i-th column of <paramref name="a"/>P
        /// was the k-th column of <paramref name="a"/>.</param>
        /// <param name="rcond"><b>rcond</b> is used to determine the effective rank of <paramref name="a"/>, which
        /// is defined as the order of the largest leading triangular
        /// submatrix R00 in the QR factorization with pivoting of <paramref name="a"/>,
        /// whose estimated condition number &lt; <b>1/rcond</b>.</param>
        /// <param name="rank"> The effective rank of <paramref name="a"/>, i.e., the order of the submatrix
        /// R00. This is the same as the order of the submatrix T00
        /// in the complete orthogonal factorization of <paramref name="a"/>.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgelsy(int m, int n, int nrhs, Complex[] a, int lda, Complex[] b, int ldb, int[] jpvt, double rcond, ref int rank);

        /// <summary>
        /// <para>Computes the minimum-norm solution to a real linear least squares problem:</para>
        /// <para>minimize || A * X - B ||</para>
        /// <para>using a complete orthogonal factorization of A. A is an M by N
        /// matrix which may be rank-deficient. Several right hand side vectors b and solution vectors x can be
        /// handled in a single call; they are stored as the columns of the M by NRHS right hand side matrix B 
        /// and the N by NRHS solution matrix X. The routine first computes a QR factorization with column pivoting:
        /// <code>
        /// A * P = Q * [ R00 R01 ]
        ///             [ 0   R11 ]
        /// </code>
        /// with R11 defined as the largest leading submatrix whose estimated
        /// condition number is less than 1/RCOND. The order of R11, RANK,
        /// is the effective rank of A.</para>
        /// <para>Then, R22 is considered to be negligible, and R12 is annihilated
        /// by orthogonal transformations from the right, arriving at the complete orthogonal factorization:
        /// <code>
        /// A * P = Q * [ T00 0 ] * Z
        ///             [ 0   0 ]
        /// </code>
        /// The minimum-norm solution is then
        /// <code>
        /// X = P * Z' [ inv(T00)*Q0'*B ]
        ///            [        0       ]
        /// </code>
        /// where Q0 consists of the first RANK columns of Q.</para>
        /// This routine is basically identical to the original xGELSX except three differences:
        /// <list>
        /// <item>The call to the subroutine xGEQPF has been substituted by the
        /// the call to the subroutine xGEQP3. This subroutine is a Blas-3
        /// version of the QR factorization with column pivoting.</item>
        /// <item>DenseMatrix B (the right hand side) is updated with Blas-3.</item>
        /// <item>The permutation of matrix B (the right hand side) is faster and more simple.</item>
        /// </list>
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of
        /// columns of matrices <paramref name="b"/> and X. <paramref name="nrhs"/> &gt;= 0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit, <b>a</b> has been overwritten by details of its complete orthogonal factorization.</param>
        /// <param name="lda"> The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="b">dimension (<paramref name="ldb"/>,<paramref name="nrhs"/>)
        /// On entry, the <paramref name="m"/> by <paramref name="nrhs"/> right hand side matrix <paramref name="b"/>.
        /// On exit, the <paramref name="n"/> by <paramref name="nrhs"/> solution matrix X.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <b>ldb</b> &gt;= max(1,<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <param name="jpvt">On entry, if <b>jpvt[i] != 0</b>, the i-th column of <paramref name="a"/> is permuted
        /// to the front of <paramref name="a"/>P, otherwise column i is a free column.
        /// On exit, if <b>jpvt[i] = k</b>, then the i-th column of <paramref name="a"/>P
        /// was the k-th column of <paramref name="a"/>.</param>
        /// <param name="rcond"><b>rcond</b> is used to determine the effective rank of <paramref name="a"/>, which
        /// is defined as the order of the largest leading triangular
        /// submatrix R00 in the QR factorization with pivoting of <paramref name="a"/>,
        /// whose estimated condition number &lt; <b>1/rcond</b>.</param>
        /// <param name="rank"> The effective rank of <paramref name="a"/>, i.e., the order of the submatrix
        /// R00. This is the same as the order of the submatrix T00
        /// in the complete orthogonal factorization of <paramref name="a"/>.</param>
        /// <param name="work">Work array. On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use.
        /// The length of <b>work</b> should be at least MN + MAX( 2*MN, NB*(<paramref name="n"/>+1), MN+MN*NB, MN+NB*<paramref name="nrhs"/> ) 
        /// where MN = min( <paramref name="m"/>, <paramref name="n"/> ) and NB is an upper bound on the blocksize returned by ILAENV
        /// for the routines ZGEQP3, ZTZRZF, ZTZRQF, ZUNMQR, and ZUNMRZ.</param>
        /// <param name="rwork">Dimension 2*<paramref name="n"/></param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgelsy(int m, int n, int nrhs, Complex[] a, int lda, Complex[] b, int ldb, int[] jpvt, double rcond, ref int rank, Complex[] work, double[] rwork);

        /// <summary>
        /// Computes a QR factorization with column pivoting of a matrix A: A*P = Q*R using Level 3 BLAS.
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit, the upper triangle of the array contains the
        /// min(<paramref name="m"/>,<paramref name="n"/>) by <paramref name="n"/> upper trapezoidal matrix R. The elements below
        /// the diagonal, together with the array <paramref name="tau"/>, represent the
        /// unitary matrix Q as a product of min(<paramref name="m"/>,<paramref name="n"/>) elementary
        /// reflectors.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="jpvt">Dimension <paramref name="n"/>.
        /// On entry, if <b>jpvt[j] != 0</b>, the j-th column of <paramref name="a"/> is permuted
        /// to the front of <paramref name="a"/>*P (a leading column); if <b>jpvt[j] == 0</b>,
        /// the j-th column of <paramref name="a"/> is a free column.
        /// On exit, if <b>jpvt[j] == k</b>, then the j-th column of <paramref name="a"/>*P was the
        /// the k-th column of <paramref name="a"/>.</param>
        /// <param name="tau">Dimension min(<paramref name="m"/>,<paramref name="n"/>). The scalar factors of the elementary reflectors.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgeqp3(int m, int n, Complex[] a, int lda, int[] jpvt, Complex[] tau);

        /// <summary>
        /// Computes a QR factorization with column pivoting of a matrix A: A*P = Q*R using Level 3 BLAS.
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit, the upper triangle of the array contains the
        /// min(<paramref name="m"/>,<paramref name="n"/>) by <paramref name="n"/> upper trapezoidal matrix R. The elements below
        /// the diagonal, together with the array <paramref name="tau"/>, represent the
        /// unitary matrix Q as a product of min(<paramref name="m"/>,<paramref name="n"/>) elementary
        /// reflectors.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="jpvt">Dimension <paramref name="n"/>.
        /// On entry, if <b>jpvt[j] != 0</b>, the j-th column of <paramref name="a"/> is permuted
        /// to the front of <paramref name="a"/>*P (a leading column); if <b>jpvt[j] == 0</b>,
        /// the j-th column of <paramref name="a"/> is a free column.
        /// On exit, if <b>jpvt[j] == k</b>, then the j-th column of <paramref name="a"/>*P was the
        /// the k-th column of <paramref name="a"/>.</param>
        /// <param name="tau">Dimension min(<paramref name="m"/>,<paramref name="n"/>). The scalar factors of the elementary reflectors.</param>
        /// <param name="work">Work array. On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use. 
        /// The length of <b>work</b> should be at least<paramref name="n"/>+1.
        /// For optimal performance the length of <b>work</b> should be at least( <paramref name="n"/>+1 )*NB, where NB is the optimal blocksize.</param>
        /// <param name="rwork">Dimension 2*<paramref name="n"/></param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgeqp3(int m, int n, Complex[] a, int lda, int[] jpvt, Complex[] tau, Complex[] work, double[] rwork);

        /// <summary>
        /// Computes an LU factorization of a general M by N matrix A
        /// using partial pivoting with row interchanges. 
        /// The factorization has the form A = P * L * U
        /// where P is a permutation matrix, L is lower triangular with unit
        /// diagonal elements (lower trapezoidal if m &gt; n), and U is upper
        /// triangular (upper trapezoidal if m &lt; n).
        /// This is the right-looking Level 3 BLAS version of the algorithm.
        /// </summary>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimensions <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/>-by-<paramref name="n"/> matrix to be factored.
        /// On exit, the factors L and U from the factorization
        /// A = P*L*U. The unit diagonal elements of L are not stored.</param>
        /// <param name="lda"> The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <param name="ipiv">Dimension min(<paramref name="m"/>,<paramref name="n"/>).
        /// The pivot indices. for 0 &lt;= i &lt;= min(<paramref name="m"/>,<paramref name="n"/>), row i of the matrix was interchanged with row <b>ipiv[i]</b>.</param>
        /// <returns>
        /// <list> 
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: U[i,i] is exactly zero. The factorization has been completed, but the factor U is exactly
        ///	singular, and division by zero will occur if it is used
        /// to solve a system of equations.</item>
        /// </list></returns>
        int Zgetrf(int m, int n, Complex[] a, int lda, out int[] ipiv);

        /// <summary>
        /// Computes the inverse of a matrix using the LU factorization
        /// computed by ZGETRF. This method inverts U and then computes 
        /// inv(A) by solving the system inv(A)*L = inv(U) for inv(A).
        /// </summary>
        /// <param name="n"> The order of the matrix <parmaref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <parmaref name="lda"/> by <parmaref name="n"/>. On entry, the factors L and U from the factorization
        /// A = P*L*U as computed by ZGETRF.
        /// On exit, the inverse of the original matrix <b>a</b>.</param>
        /// <param name="lda">The leading dimension of the array <parmaref name="a"/>. <b>lda</b> &gt;= max(1,<parmaref name="n"/>).</param>
        /// <param name="ipiv">Dimension <parmaref name="n"/>. The pivot indices from ZGETRF; for 0 &lt;= i &lt;= <parmaref name="n"/>-1, row i of the
        /// matrix was interchanged with row <b>ipiv[i]</b>.</param>
        /// <returns><list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value</item>
        /// <item>&gt; 0: U[i,i] is exactly zero; the matrix is
        /// singular and its inverse could not be computed.</item>
        /// </list></returns>
        int Zgetri(int n, Complex[] a, int lda, int[] ipiv);

        /// <summary>
        /// Computes the inverse of a matrix using the LU factorization
        /// computed by ZGETRF. This method inverts U and then computes 
        /// inv(A) by solving the system inv(A)*L = inv(U) for inv(A).
        /// </summary>
        /// <param name="n"> The order of the matrix <parmaref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <parmaref name="lda"/> by <parmaref name="n"/>. On entry, the factors L and U from the factorization
        /// A = P*L*U as computed by ZGETRF.
        /// On exit, the inverse of the original matrix <b>a</b>.</param>
        /// <param name="lda">The leading dimension of the array <parmaref name="a"/>. <b>lda</b> &gt;= max(1,<parmaref name="n"/>).</param>
        /// <param name="ipiv">Dimension <parmaref name="n"/>. The pivot indices from ZGETRF; for 0 &lt;= i &lt;= <parmaref name="n"/>-1, row i of the
        /// matrix was interchanged with row <b>ipiv[i]</b>.</param>
        /// <param name="work">Work array. The length of <b>work</b> should be at least <paramref name="n"/>*NB,
        /// where NB is the optimal blocksize returned by ILAENV.
        /// On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use.</param>
        /// <returns><list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value</item>
        /// <item>&gt; 0: U[i,i] is exactly zero; the matrix is
        /// singular and its inverse could not be computed.</item>
        /// </list></returns>
        int Zgetri(int n, Complex[] a, int lda, int[] ipiv, Complex[] work);

        /// <summary>
        /// Solves a system of linear equations <code>A * X = B, A**T * X = B, or A**H * X = B</code>
        /// with a general N by N matrix A using the LU factorization computed by ZGETRF.
        /// </summary>
        /// <param name="trans"> Specifies the form of the system of equations:
        /// <list>
        /// <item>'N': A * X = B (No transpose)</item>
        /// <item>'T': A**T * X = B (Transpose)</item>
        /// <item>'C': A**H * X = B (Conjugate transpose)</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns
        /// of the matrix <paramref name="b"/>. <b>nrhs &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// The factors L and U from the factorization A = P*L*U as computed by ZGETRF.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <param name="ipiv">Dimension <parmaref name="n"/>. The pivot indices from ZGETRF; for 0 &lt;= i &lt;= <parmaref name="n"/>-1, row i of the
        /// matrix was interchanged with row <b>ipiv[i]</b>.</param>
        /// <param name="b">Dimension <parmaref name="ldb"/> by <parmaref name="nrhs"/>.
        /// On entry, the right hand side matrix <b>b</b>.
        /// On exit, the solution matrix X.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <b>ldb</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgetrs(char trans, int n, int nrhs, Complex[] a, int lda, int[] ipiv, Complex[] b, int ldb);

        /// <summary>
        /// Computes the Cholesky factorization of a Hermitian
        /// positive definite matrix A.
        /// The factorization has the form 
        /// <code>
        /// A = U**H * U, if UPLO == 'U' 
        /// A = L * L**H, if UPLO == 'L'
        /// </code>
        /// where U is an upper triangular matrix and L is lower triangular.
        /// This is the block version of the algorithm, calling Level 3 BLAS.
        /// </summary>
        /// <param name="uplo"><list>
        /// <item>'U': Upper triangle of A is stored.</item>
        /// <item>'L': Lower triangle of A is stored.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// On entry, the Hermitian matrix <b>a</b>. 
        /// If <paramref name="uplo"/> == 'U', the leading <paramref name="n"/> by <paramref name="n"/> upper triangular part of <b>a</b>
        /// contains the upper triangular part of the matrix <b>a</b>, and the strictly lower triangular part of <b>a</b> is not referenced.
        /// If <paramref name="uplo"/> == 'L', the leading <paramref name="n"/> by <paramref name="n"/> lower triangular part of <b>a</b>
        /// contains the lower triangular part of the matrix <b>a</b>, and the strictly upper triangular part of <b>a</b> is not referenced.
        /// On exit, the factor U or L from the Cholesky factorization <b>A = U**H*U or A = L*L**H</b>.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <returns><list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: the leading minor of order i is not positive definite, and the factorization could not be completed.</item>
        /// </list></returns>
        int Zpotrf(char uplo, int n, Complex[] a, int lda);

        /// <summary>
        /// Computes the inverse of a Hermitian positive definite
        /// matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
        /// computed by ZPOTRF.
        /// </summary>
        /// <param name="uplo"><list>
        /// <item>'U': Upper triangle of A is stored.</item>
        /// <item>'L': Lower triangle of A is stored.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// On entry, the triangular factor U or L from the Cholesky
        /// factorization A = U**H*U or A = L*L**H, as computed by ZPOTRF.
        /// On exit, the upper or lower triangle of the (Hermitian) inverse of <b>a</b>, 
        /// overwriting the input factor U or L.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <returns><list>
        /// <item>0: successful exit</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: the [i,i] element of the factor U or L is zero, and the inverse could not be computed.</item>
        /// </list></returns>
        int Zpotri(char uplo, int n, Complex[] a, int lda);

        /// <summary>
        /// Solves a system of linear equations A*X = B with a Hermitian
        /// positive definite matrix A using the Cholesky factorization
        /// <code>A = U**H*U or A = L*L**H</code> computed by ZPOTRF.
        /// </summary>
        /// <param name="uplo"><list>
        /// <item>'U': Upper triangle of A is stored.</item>
        /// <item>'L': Lower triangle of A is stored.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="nrhs">The number of right hand sides, i.e., the number of columns
        /// of the matrix <paramref name="b"/>. <b>nrhs &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// The triangular factor U or L from the Cholesky factorization
        /// A = U**H*U or A = L*L**H, as computed by ZPOTRF.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <param name="b">Dimension <paramref name="ldb"/> by <paramref name="nrhs"/>.
        /// On entry, the right hand side matrix <b>b</b>.
        /// On exit, the solution matrix X.</param>
        /// <param name="ldb">The leading dimension of the array <paramref name="b"/>. <b>ldb</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zpotrs(char uplo, int n, int nrhs, Complex[] a, int lda, Complex[] b, int ldb);

        /// <summary>
        /// Generates an M by N real matrix Q with orthonormal columns,
        /// which is defined as the first N columns of a product of K elementary
        /// reflectors of order M: <code>Q = H(0) H(1) . . . H(k-1)</code>
        /// as returned by ZGEQRF.
        /// </summary>
        /// <param name="m">The number of rows of the matrix Q. <b>m &gt;= 0</b>.</param>
        /// <param name="n">The number of columns of the matrix Q. <b>m &gt;= n &gt;= 0</b></param>
        /// <param name="k">The number of elementary reflectors whose product defines the matrix Q. <b>n &gt;= k &gt;= 0</b>.</param>
        /// <param name="a">Dimension LDA by N.
        /// On entry, the i-th column must contain the vector which defines the elementary reflector H(0), for i = 0,2,...,k-1, as
        /// returned by ZGEQRF in the first k columns of its array argument <b>a</b>.
        /// On exit, the <parmaref name="m"/> by <parmaref name="n"/> matrix Q.</param>
        /// <param name="lda">The first dimension of the array <parmaref name="a"/> . <b>lda</b> &gt;= max(1,<parmaref name="m"/>).</param>
        /// <param name="tau">Dimension K. <b>tau[i]</b> must contain the scalar factor of the elementary
        /// reflector H(i), as returned by ZGEQRF.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zungqr(int m, int n, int k, Complex[] a, int lda, Complex[] tau);

        /// <summary>
        /// Generates an M by N real matrix Q with orthonormal columns,
        /// which is defined as the first N columns of a product of K elementary
        /// reflectors of order M: <code>Q = H(0) H(1) . . . H(k-1)</code>
        /// as returned by ZGEQRF.
        /// </summary>
        /// <param name="m">The number of rows of the matrix Q. <b>m &gt;= 0</b>.</param>
        /// <param name="n">The number of columns of the matrix Q. <b>m &gt;= n &gt;= 0</b></param>
        /// <param name="k">The number of elementary reflectors whose product defines the matrix Q. <b>n &gt;= k &gt;= 0</b>.</param>
        /// <param name="a">Dimension LDA by N.
        /// On entry, the i-th column must contain the vector which defines the elementary reflector H(0), for i = 0,2,...,k-1, as
        /// returned by ZGEQRF in the first k columns of its array argument <b>a</b>.
        /// On exit, the <parmaref name="m"/> by <parmaref name="n"/> matrix Q.</param>
        /// <param name="lda">The first dimension of the array <parmaref name="a"/> . <b>lda</b> &gt;= max(1,<parmaref name="m"/>).</param>
        /// <param name="tau">Dimension K. <b>tau[i]</b> must contain the scalar factor of the elementary
        /// reflector H(i), as returned by ZGEQRF.</param>
        /// <param name="work">Work array. The length of <b>work</b> should be at least <paramref name="n"/>*NB,
        /// where NB is the optimal blocksize returned by ILAENV.
        /// On exit, <b>work[0]</b> returns the optimal length of <b>work</b> to use.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zungqr(int m, int n, int k, Complex[] a, int lda, Complex[] tau, Complex[] work);

        /// <summary>
        /// <para>Overwrites the general real M-by-N matrix C with
        /// <pre>             SIDE = 'L'    SIDE = 'R'</pre>
        /// <pre>TRANS = 'N':   Q * C         C * Q</pre>
        /// <pre>TRANS = 'T':   Q**H * C      C * Q**H</pre>																				
        /// where Q is a real orthogonal matrix defined as the product of k elementary reflectors
        /// Q = H(0) H(1) . . . H(k-1) as returned by ZGEQRF. Q is of order M if 
        /// SIDE == 'L' and of order N if SIDE == 'R'.</para>
        /// </summary>
        /// <param name="side"><list>
        /// <item>'L': apply Q or Q**H from the Left.</item>
        /// <item>'R': apply Q or Q**H from the Right.</item>
        /// </list></param>
        /// <param name="trans"><list>
        /// <item>'N': No transpose, apply Q;</item>
        /// <item>'T': Transpose, apply Q**H.</item>
        /// </list></param>
        /// <param name="m">The number of rows of the matrix <paramref name="c"/>. <b>m &gt;= 0</b>.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="c"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="k"> The number of elementary reflectors whose product defines the matrix Q.
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <paramref name="m"/> &gt;= <paramref name="k"/> &gt;= 0.</item>
        /// <item>If <paramref name="side"/> == 'R', <paramref name="n"/> &gt;= <paramref name="k"/> &gt;= 0.</item>
        /// </list></param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="k"/>.
        /// The i-th column must contain the vector which defines the
        /// elementary reflector H(i), for i = 0,2,...,<paramref name="k"/>-1, as returned by
        ///	ZGEQRF in the first <paramref name="k"/> columns of its array argument <paramref name="a"/>.
        ///	<paramref name="a"/> is modified by the routine but restored on exit.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/> .
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <b>lda</b> &gt;= max(1,<paramref name="m"/>).</item>
        /// <item>If <paramref name="side"/> == 'R', <b>lda</b> &gt;= max(1,<paramref name="n"/>).</item>
        /// </list></param>
        /// <param name="tau">Dimension <paramref name="k"/>. <b>tau[i]</b> must contain the scalar factor of the elementary
        /// reflector H(i), as returned by ZGEQRF.</param>
        /// <param name="c">dimension <paramref name="ldc"/> by <paramref name="n"/>. 
        /// On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>c</b>.
        /// On exit, <b>c</b> is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.</param>
        /// <param name="ldc">The leading dimension of the array <paramref name="c"/>.<b>ldc</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zunmqr(char side, char trans, int m, int n, int k, Complex[] a, int lda, Complex[] tau, Complex[] c, int ldc);

        /// <summary>
        /// <para>Overwrites the general real M-by-N matrix C with
        /// <pre>             SIDE = 'L'    SIDE = 'R'</pre>
        /// <pre>TRANS = 'N':   Q * C         C * Q</pre>
        /// <pre>TRANS = 'T':   Q**H * C      C * Q**H</pre>																			
        /// where Q is a real orthogonal matrix defined as the product of k elementary reflectors
        /// Q = H(0) H(1) . . . H(k-1) as returned by ZGEQRF. Q is of order M if 
        /// SIDE == 'L' and of order N if SIDE == 'R'.</para>
        /// </summary>
        /// <param name="side"><list>
        /// <item>'L': apply Q or Q**H from the Left.</item>
        /// <item>'R': apply Q or Q**H from the Right.</item>
        /// </list></param>
        /// <param name="trans"><list>
        /// <item>'N': No transpose, apply Q;</item>
        /// <item>'T': Transpose, apply Q**H.</item>
        /// </list></param>
        /// <param name="m">The number of rows of the matrix <paramref name="c"/>. <b>m &gt;= 0</b>.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="c"/>. <b>n &gt;= 0</b>.</param>
        /// <param name="k"> The number of elementary reflectors whose product defines the matrix Q.
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <paramref name="m"/> &gt;= <paramref name="k"/> &gt;= 0.</item>
        /// <item>If <paramref name="side"/> == 'R', <paramref name="n"/> &gt;= <paramref name="k"/> &gt;= 0.</item>
        /// </list></param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="k"/>.
        /// The i-th column must contain the vector which defines the
        /// elementary reflector H(i), for i = 0,2,...,<paramref name="k"/>-1, as returned by
        ///	ZGEQRF in the first <paramref name="k"/> columns of its array argument <paramref name="a"/>.
        ///	<paramref name="a"/> is modified by the routine but restored on exit.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/> .
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <b>lda</b> &gt;= max(1,<paramref name="m"/>).</item>
        /// <item>If <paramref name="side"/> == 'R', <b>lda</b> &gt;= max(1,<paramref name="n"/>).</item>
        /// </list></param>
        /// <param name="tau">Dimension <paramref name="k"/>. <b>tau[i]</b> must contain the scalar factor of the elementary
        /// reflector H(i), as returned by ZGEQRF.</param>
        /// <param name="c">dimension <paramref name="ldc"/> by <paramref name="n"/>. 
        /// On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>c</b>.
        /// On exit, <b>c</b> is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.</param>
        /// <param name="ldc">The leading dimension of the array <paramref name="c"/>.<b>ldc</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="work">Work array. The length of <b>work</b> should be:
        /// <list>
        /// <item>If <paramref name="side"/> == 'L', <paramref name="n"/> * NB.</item>
        /// <item>if <paramref name="side"/> == 'R', <paramref name="m"/> * NB.</item>
        /// </list>
        /// where NB is the optimal blocksize returned by ILAENV.
        /// On exit, <b>work[0]</b> returns the optimal length of <b>work</b> to use.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zunmqr(char side, char trans, int m, int n, int k, Complex[] a, int lda, Complex[] tau, Complex[] c, int ldc, Complex[] work);

        /// <summary>
        /// Computes the singular value decomposition (SVD) of a complex
        /// M by N matrix A, optionally computing the left and/or right singular
        /// vectors. The SVD is written A = U * SIGMA * conjugate-transpose(V)
        /// where SIGMA is an M by N matrix which is zero except for its
        /// min(m,n) diagonal elements, U is an M by M unitary matrix, and V is an N-by-N unitary matrix. The diagonal elements of SIGMA
        /// are the singular values of A. They are real and non-negative, and
        /// are returned in descending order. The first min(m,n) columns of
        /// U and V are the left and right singular vectors of A.
        /// Note that the routine returns V**H, not V.
        /// </summary>
        /// <param name="jobu">
        /// Specifies options for computing all or part of the matrix <paramref name="u"/>:
        /// <list>
        /// <item> = '<paramref name="a"/>': all <paramref name="m"/> columns of <paramref name="u"/> are returned in array <paramref name="u"/>.</item>
        /// <item> = '<paramref name="s"/>': the first min(m,n) columns of <paramref name="u"/> (the left singular vectors) are returned in the array <paramref name="u"/>.</item>
        /// <item> = 'O': the first min(m,n) columns of <paramref name="u"/> (the left singular vectors) are overwritten on the array <paramref name="a"/>.</item>
        /// <item> = '<paramref name="n"/>': no columns of <paramref name="u"/> (no left singular vectors) are computed. </item>
        /// </list>
        /// </param>
        /// <param name="jobvt">
        /// Specifies options for computing all or part of the matrix V**T:
        /// <list>
        /// <item>'<paramref name="a"/>': all <paramref name="n"/> rows of V**T are returned in the array <paramref name="vt"/>.</item>
        /// <item>'<paramref name="s"/>': the first min(m,n) rows of V**T (the right singular vectors) are returned in the array <paramref name="vt"/>.</item>
        /// <item>'O': the first min(m,n) rows of V**T (the right singular vectors) are overwritten on the array <paramref name="a"/>.</item>
        /// <item>'<paramref name="n"/>': no rows of V**T (no right singular vectors) are computed.</item>
        /// </list>
        /// </param>
        /// <param name="m">The number of rows of the input matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the input matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit:<list>
        /// <item>if <paramref name="jobu"/> == 'O', <b>a</b> is overwritten with the first min(m,n) columns of <paramref name="u"/> (the left singular vectors, stored column wise).</item>
        /// <item>if <paramref name="jobvt"/> == 'O', <b>a</b> is overwritten with the first min(m,n)rows of V**T (the right singular vectors, stored row wise).</item>
        /// <item>if <paramref name="jobu"/> != 'O' and <paramref name="jobvt"/> != 'O', the contents of <b>a</b> are destroyed.</item>
        /// </list></param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <param name="s">Dimension min(<paramref name="m"/>,<paramref name="n"/>). The singular values of <paramref name="a"/>, sorted so that <b>s[i] &gt;= s[i+1]</b>.</param>
        /// <param name="u">Dimension <paramref name="ldu"/> by UCOL, <paramref name="ldu"/> by <paramref name="m"/> if <paramref name="jobu"/> = '<paramref name="a"/>', or <paramref name="ldu"/> by min(<paramref name="m"/>,<paramref name="n"/>) if <paramref name="jobu"/> = '<paramref name="s"/>'.
        /// <list>
        /// <item>If <paramref name="jobu"/> == '<paramref name="a"/>', <b>u</b> contains the <paramref name="m"/> by <paramref name="m"/> orthogonal matrix <b>u</b>.</item>
        /// <item>If <paramref name="jobu"/> == '<paramref name="s"/>', <b>u</b> contains the first min(m,n) columns of <b>u</b> (the left singular vectors, stored column wise).</item>
        /// <item>if <paramref name="jobu"/> == '<paramref name="n"/>' or 'O', <b>u</b> is not referenced.</item>
        /// </list></param>
        /// <param name="ldu">The leading dimension of the array <paramref name="u"/>. <b>ldu &gt;= 1</b>. If <paramref name="jobu"/> = '<paramref name="s"/>' or '<paramref name="a"/>', <b>ldu &gt;= <paramref name="m"/></b>.</param>
        /// <param name="vt">Dimension <paramref name="ldvt"/> by <paramref name="n"/>.
        /// <list>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="a"/>', <b>vt</b> contains the <paramref name="n"/> by <paramref name="n"/> orthogonal matrix V**T.</item>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="s"/>', <b>vt</b> contains the first min(m,n) rows of V**T (the right singular vectors, stored row wise).</item>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="n"/>' or 'O', <b>vt</b> is not referenced.</item>
        /// </list></param>
        /// <param name="ldvt">The leading dimension of the array <paramref name="vt"/>. <b>ldvt &gt;= 1</b>. If <paramref name="jobvt"/> == '<paramref name="a"/>', 
        /// <b>ldvt &gt;= <paramref name="n"/></b>; if <paramref name="jobvt"/> == '<paramref name="s"/>', <b><paramref name="ldvt"/> &gt;= min(<paramref name="m"/>,<paramref name="n"/>)</b>.</param>
        /// <returns><list>
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: if ZBDSQR did not converge, specifies how many
        /// superdiagonals of an intermediate bidiagonal form B
        /// did not converge to zero.</item>
        /// </list></returns>
        /// <remarks><paramref name="jobvt"/> and <paramref name="jobu"/> cannot both be 'O'.</remarks>
        int Zgesvd(char jobu, char jobvt, int m, int n, Complex[] a, int lda, double[] s, Complex[] u, int ldu, Complex[] vt, int ldvt);

        /// <summary>
        /// Computes the singular value decomposition (SVD) of a complex
        /// M by N matrix A, optionally computing the left and/or right singular
        /// vectors. The SVD is written A = U * SIGMA * conjugate-transpose(V)
        /// where SIGMA is an M by N matrix which is zero except for its
        /// min(m,n) diagonal elements, U is an M by M unitary matrix, and V is an N-by-N unitary matrix. The diagonal elements of SIGMA
        /// are the singular values of A. They are real and non-negative, and
        /// are returned in descending order. The first min(m,n) columns of
        /// U and V are the left and right singular vectors of A.
        /// Note that the routine returns V**H, not V.
        /// </summary>
        /// <param name="jobu">
        /// Specifies options for computing all or part of the matrix <paramref name="u"/>:
        /// <list>
        /// <item> = '<paramref name="a"/>': all <paramref name="m"/> columns of <paramref name="u"/> are returned in array <paramref name="u"/>.</item>
        /// <item> = '<paramref name="s"/>': the first min(m,n) columns of <paramref name="u"/> (the left singular vectors) are returned in the array <paramref name="u"/>.</item>
        /// <item> = 'O': the first min(m,n) columns of <paramref name="u"/> (the left singular vectors) are overwritten on the array <paramref name="a"/>.</item>
        /// <item> = '<paramref name="n"/>': no columns of <paramref name="u"/> (no left singular vectors) are computed. </item>
        /// </list>
        /// </param>
        /// <param name="jobvt">
        /// Specifies options for computing all or part of the matrix V**T:
        /// <list>
        /// <item>'<paramref name="a"/>': all <paramref name="n"/> rows of V**T are returned in the array <paramref name="vt"/>.</item>
        /// <item>'<paramref name="s"/>': the first min(m,n) rows of V**T (the right singular vectors) are returned in the array <paramref name="vt"/>.</item>
        /// <item>'O': the first min(m,n) rows of V**T (the right singular vectors) are overwritten on the array <paramref name="a"/>.</item>
        /// <item>'<paramref name="n"/>': no rows of V**T (no right singular vectors) are computed.</item>
        /// </list>
        /// </param>
        /// <param name="m">The number of rows of the input matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the input matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. On entry, the <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.
        /// On exit:<list>
        /// <item>if <paramref name="jobu"/> == 'O', <b>a</b> is overwritten with the first min(m,n) columns of <paramref name="u"/> (the left singular vectors, stored column wise).</item>
        /// <item>if <paramref name="jobvt"/> == 'O', <b>a</b> is overwritten with the first min(m,n)rows of V**T (the right singular vectors, stored row wise).</item>
        /// <item>if <paramref name="jobu"/> != 'O' and <paramref name="jobvt"/> != 'O', the contents of <b>a</b> are destroyed.</item>
        /// </list></param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <param name="s">Dimension min(<paramref name="m"/>,<paramref name="n"/>). The singular values of <paramref name="a"/>, sorted so that <b>s[i] &gt;= s[i+1]</b>.</param>
        /// <param name="u">Dimension <paramref name="ldu"/> by UCOL, <paramref name="ldu"/> by <paramref name="m"/> if <paramref name="jobu"/> = '<paramref name="a"/>', or <paramref name="ldu"/> by min(<paramref name="m"/>,<paramref name="n"/>) if <paramref name="jobu"/> = '<paramref name="s"/>'.
        /// <list>
        /// <item>If <paramref name="jobu"/> == '<paramref name="a"/>', <b>u</b> contains the <paramref name="m"/> by <paramref name="m"/> orthogonal matrix <b>u</b>.</item>
        /// <item>If <paramref name="jobu"/> == '<paramref name="s"/>', <b>u</b> contains the first min(m,n) columns of <b>u</b> (the left singular vectors, stored column wise).</item>
        /// <item>if <paramref name="jobu"/> == '<paramref name="n"/>' or 'O', <b>u</b> is not referenced.</item>
        /// </list></param>
        /// <param name="ldu">The leading dimension of the array <paramref name="u"/>. <b>ldu &gt;= 1</b>. If <paramref name="jobu"/> = '<paramref name="s"/>' or '<paramref name="a"/>', <b>ldu &gt;= <paramref name="m"/></b>.</param>
        /// <param name="vt">Dimension <paramref name="ldvt"/> by <paramref name="n"/>.
        /// <list>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="a"/>', <b>vt</b> contains the <paramref name="n"/> by <paramref name="n"/> orthogonal matrix V**T.</item>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="s"/>', <b>vt</b> contains the first min(m,n) rows of V**T (the right singular vectors, stored row wise).</item>
        /// <item>If <paramref name="jobvt"/> == '<paramref name="n"/>' or 'O', <b>vt</b> is not referenced.</item>
        /// </list></param>
        /// <param name="ldvt">The leading dimension of the array <paramref name="vt"/>. <b>ldvt &gt;= 1</b>. If <paramref name="jobvt"/> == '<paramref name="a"/>', 
        /// <b>ldvt &gt;= <paramref name="n"/></b>; if <paramref name="jobvt"/> == '<paramref name="s"/>', <b><paramref name="ldvt"/> &gt;= min(<paramref name="m"/>,<paramref name="n"/>)</b>.</param>
        /// <param name="work">Work array. On exit, <paramref name="work"/>[0] returns the optimal length of <paramref name="work"/> to use.
        /// The length of <b>work</b> should be at least 2*min(<paramref name="m"/>,<paramref name="n"/>)+max(<paramref name="m"/>,<paramref name="n"/>).</param>
        /// <param name="rwork">Dimension 5*min(<paramref name="m"/>,<paramref name="n"/>).
        /// On exit, if method returns a value greater than 0, <b>rwork[0:min(<paramref name="m"/>,<paramref name="n"/>)-2]</b> contains the unconverged superdiagonal elements of an upper bidiagonal
        /// matrix B whose diagonal is in <paramref name="s"/> (not necessarily sorted).
        /// B satisfies <paramref name="a"/> = <paramref name="u"/> * B * <paramref name="vt"/>, so it has the same singular
        /// values as <paramref name="a"/>, and singular vectors related by <paramref name="u"/> and <paramref name="vt"/>.</param>
        /// <returns><list>
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// <item>&gt; 0: if ZBDSQR did not converge, specifies how many
        /// superdiagonals of an intermediate bidiagonal form B
        /// did not converge to zero.See the description of <paramref name="rwork"/>
        /// above for details.</item>
        /// </list></returns>
        /// <remarks><paramref name="jobvt"/> and <paramref name="jobu"/> cannot both be 'O'.</remarks>
        int Zgesvd(char jobu, char jobvt, int m, int n, Complex[] a, int lda, double[] s, Complex[] u, int ldu, Complex[] vt, int ldvt, Complex[] work, double[] rwork);

        /// <summary>
        /// Returns the value of the one norm, or the Frobenius norm, or
        /// the infinity norm, or the element of largest absolute value of a
        /// matrix A.
        /// </summary>
        /// <param name="norm">Specifies the value to be returned.</param>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <b>m &gt;= 0</b>. When <b>m == 0</b>,
        /// 0 is returned.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>. When <b>n == 0</b>,
        /// 0 is returned.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. The <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <returns><list>
        /// <item>max(abs(<paramref name="a"/>[i,j])), if <paramref name="norm"/> == 'M' or 'm'.</item>
        /// <item>one norm of a matrix (maximum column sum), if <paramref name="norm"/> == '1', 'O' or 'o'.</item>
        /// <item>infinity norm of a matrix (maximum row sum), if <paramref name="norm"/> == 'I' or 'i'.</item>
        /// <item>Frobenius norm of a matrix (square root of sum of squares), if <paramref name="norm"/> == 'F', 'f', 'E' or 'e'.</item>
        /// </list>
        /// Note that max(abs(<paramref name="a"/>[i,j])) is not a matrix norm. 
        ///</returns>
        double Dlange(char norm, int m, int n, double[] a, int lda);

        /// <summary>
        /// Returns the value of the one norm, or the Frobenius norm, or
        /// the infinity norm, or the element of largest absolute value of a
        /// matrix A.
        /// </summary>
        /// <param name="norm">Specifies the value to be returned.</param>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <b>m &gt;= 0</b>. When <b>m == 0</b>,
        /// 0 is returned.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>. When <b>n == 0</b>,
        /// 0 is returned.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. The <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <param name="work">Work array with a length of at least <paramref name="m"/> when <paramref name="norm"/> == 'I',
        /// otherwise, the array is not referenced.</param>
        /// <returns><list>
        /// <item>max(abs(<paramref name="a"/>[i,j])), if <paramref name="norm"/> == 'M' or 'm'.</item>
        /// <item>one norm of a matrix (maximum column sum), if <paramref name="norm"/> == '1', 'O' or 'o'.</item>
        /// <item>infinity norm of a matrix (maximum row sum), if <paramref name="norm"/> == 'I' or 'i'.</item>
        /// <item>Frobenius norm of a matrix (square root of sum of squares), if <paramref name="norm"/> == 'F', 'f', 'E' or 'e'.</item>
        /// </list>
        /// Note that max(abs(<paramref name="a"/>[i,j])) is not a matrix norm. 
        ///</returns>
        double Dlange(char norm, int m, int n, double[] a, int lda, double[] work);

        /// <summary>
        /// Returns the value of the one norm, or the Frobenius norm, or
        /// the infinity norm, or the element of largest absolute value of a
        /// matrix A.
        /// </summary>
        /// <param name="norm">Specifies the value to be returned.</param>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <b>m &gt;= 0</b>. When <b>m == 0</b>,
        /// 0 is returned.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>. When <b>n == 0</b>,
        /// 0 is returned.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. The <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <returns><list>
        /// <item>max(abs(<paramref name="a"/>[i,j])), if <paramref name="norm"/> == 'M' or 'm'.</item>
        /// <item>one norm of a matrix (maximum column sum), if <paramref name="norm"/> == '1', 'O' or 'o'.</item>
        /// <item>infinity norm of a matrix (maximum row sum), if <paramref name="norm"/> == 'I' or 'i'.</item>
        /// <item>Frobenius norm of a matrix (square root of sum of squares), if <paramref name="norm"/> == 'F', 'f', 'E' or 'e'.</item>
        /// </list>
        /// Note that max(abs(<paramref name="a"/>[i,j])) is not a matrix norm. 
        ///</returns>
        double Zlange(char norm, int m, int n, Complex[] a, int lda);

        /// <summary>
        /// Returns the value of the one norm, or the Frobenius norm, or
        /// the infinity norm, or the element of largest absolute value of a
        /// matrix A.
        /// </summary>
        /// <param name="norm">Specifies the value to be returned.</param>
        /// <param name="m">The number of rows of the matrix <paramref name="a"/>. <b>m &gt;= 0</b>. When <b>m == 0</b>,
        /// 0 is returned.</param>
        /// <param name="n">The number of columns of the matrix <paramref name="a"/>. <b>n &gt;= 0</b>. When <b>n == 0</b>,
        /// 0 is returned.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>. The <paramref name="m"/> by <paramref name="n"/> matrix <b>a</b>.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>. <b>lda &gt;= max(1,<paramref name="m"/>)</b>.</param>
        /// <param name="work">Work array with a length of at least <paramref name="m"/> when <paramref name="norm"/> == 'I',
        /// otherwise, the array is not referenced.</param>
        /// <returns><list>
        /// <item>max(abs(<paramref name="a"/>[i,j])), if <paramref name="norm"/> == 'M' or 'm'.</item>
        /// <item>one norm of a matrix (maximum column sum), if <paramref name="norm"/> == '1', 'O' or 'o'.</item>
        /// <item>infinity norm of a matrix (maximum row sum), if <paramref name="norm"/> == 'I' or 'i'.</item>
        /// <item>Frobenius norm of a matrix (square root of sum of squares), if <paramref name="norm"/> == 'F', 'f', 'E' or 'e'.</item>
        /// </list>
        /// Note that max(abs(<paramref name="a"/>[i,j])) is not a matrix norm. 
        ///</returns>
        double Zlange(char norm, int m, int n, Complex[] a, int lda, double[] work);

        /// <summary>
        ///  Estimates the reciprocal of the condition number of a general
        ///  real matrix A, in either the 1-norm or the infinity-norm, using
        ///  the LU factorization computed by DGETRF.
        ///  An estimate is obtained for norm(inv(A)), and the reciprocal of the
        ///  condition number is computed as RCOND = 1 / ( norm(A) * norm(inv(A)) ).
        /// </summary>
        /// <param name="norm">Specifies whether the 1-norm condition number or the infinity-norm condition number is required:
        /// <list>
        /// <item>'1' or 'O':  1-norm;</item>
        /// <item>'I': Infinity-norm.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>.  <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// The factors L and U from the factorization A = P*L*U as computed by DGETRF.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>.  <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <param name="anorm"><list>
        /// <item>If <paramref name="norm"/> == '1' or 'O', the 1-norm of the original matrix A.</item>
        /// <item>If <paramref name="norm"/> == 'I', the infinity-norm of the original matrix A.</item>
        /// </list></param>
        /// <param name="rcond">The reciprocal of the condition number of the matrix <paramref name="a"/>,
        /// computed as 1/(norm(<paramref name="a"/>) * norm(inv(<paramref name="a"/>))).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgecon(char norm, int n, double[] a, int lda, double anorm, ref double rcond);

        /// <summary>
        ///  Estimates the reciprocal of the condition number of a general
        ///  real matrix A, in either the 1-norm or the infinity-norm, using
        ///  the LU factorization computed by DGETRF.
        ///  An estimate is obtained for norm(inv(A)), and the reciprocal of the
        ///  condition number is computed as RCOND = 1 / ( norm(A) * norm(inv(A)) ).
        /// </summary>
        /// <param name="norm">Specifies whether the 1-norm condition number or the infinity-norm condition number is required:
        /// <list>
        /// <item>'1' or 'O':  1-norm;</item>
        /// <item>'I': Infinity-norm.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>.  <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// The factors L and U from the factorization A = P*L*U as computed by DGETRF.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>.  <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <param name="anorm"><list>
        /// <item>If <paramref name="norm"/> == '1' or 'O', the 1-norm of the original matrix A.</item>
        /// <item>If <paramref name="norm"/> == 'I', the infinity-norm of the original matrix A.</item>
        /// </list></param>
        /// <param name="rcond">The reciprocal of the condition number of the matrix <paramref name="a"/>,
        /// computed as 1/(norm(<paramref name="a"/>) * norm(inv(<paramref name="a"/>))).</param>
        /// <param name="work">Dimension 4*<paramref name="n"/>.</param>
        /// <param name="iwork">Dimension <paramref name="n"/>.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Dgecon(char norm, int n, double[] a, int lda, double anorm, ref double rcond, double[] work, int[] iwork);

        /// <summary>
        ///  Estimates the reciprocal of the condition number of a general
        ///  real matrix A, in either the 1-norm or the infinity-norm, using
        ///  the LU factorization computed by ZGETRF.
        ///  An estimate is obtained for norm(inv(A)), and the reciprocal of the
        ///  condition number is computed as RCOND = 1 / ( norm(A) * norm(inv(A)) ).
        /// </summary>
        /// <param name="norm">Specifies whether the 1-norm condition number or the infinity-norm condition number is required:
        /// <list>
        /// <item>'1' or 'O':  1-norm;</item>
        /// <item>'I': Infinity-norm.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>.  <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// The factors L and U from the factorization A = P*L*U as computed by ZGETRF.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>.  <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <param name="anorm"><list>
        /// <item>If <paramref name="norm"/> == '1' or 'O', the 1-norm of the original matrix A.</item>
        /// <item>If <paramref name="norm"/> == 'I', the infinity-norm of the original matrix A.</item>
        /// </list></param>
        /// <param name="rcond">The reciprocal of the condition number of the matrix <paramref name="a"/>,
        /// computed as 1/(norm(<paramref name="a"/>) * norm(inv(<paramref name="a"/>))).</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgecon(char norm, int n, Complex[] a, int lda, double anorm, ref double rcond);

        /// <summary>
        ///  Estimates the reciprocal of the condition number of a general
        ///  real matrix A, in either the 1-norm or the infinity-norm, using
        ///  the LU factorization computed by ZGETRF.
        ///  An estimate is obtained for norm(inv(A)), and the reciprocal of the
        ///  condition number is computed as RCOND = 1 / ( norm(A) * norm(inv(A)) ).
        /// </summary>
        /// <param name="norm">Specifies whether the 1-norm condition number or the infinity-norm condition number is required:
        /// <list>
        /// <item>'1' or 'O':  1-norm;</item>
        /// <item>'I': Infinity-norm.</item>
        /// </list></param>
        /// <param name="n">The order of the matrix <paramref name="a"/>.  <b>n &gt;= 0</b>.</param>
        /// <param name="a">Dimension <paramref name="lda"/> by <paramref name="n"/>.
        /// The factors L and U from the factorization A = P*L*U as computed by ZGETRF.</param>
        /// <param name="lda">The leading dimension of the array <paramref name="a"/>.  <b>lda</b> &gt;= max(1,<paramref name="n"/>).</param>
        /// <param name="anorm"><list>
        /// <item>If <paramref name="norm"/> == '1' or 'O', the 1-norm of the original matrix A.</item>
        /// <item>If <paramref name="norm"/> == 'I', the infinity-norm of the original matrix A.</item>
        /// </list></param>
        /// <param name="rcond">The reciprocal of the condition number of the matrix <paramref name="a"/>,
        /// computed as 1/(norm(<paramref name="a"/>) * norm(inv(<paramref name="a"/>))).</param>
        /// <param name="work">Dimension 2*<paramref name="n"/>.</param>
        /// <param name="rwork">Dimension 2*<paramref name="n"/>.</param>
        /// <returns>0 on successful exit, or -i, the i-th argument had an illegal value.</returns>
        int Zgecon(char norm, int n, Complex[] a, int lda, double anorm, ref double rcond, Complex[] work, double[] rwork);

        /// <summary>
        /// ILAENV is called from the LAPACK routines to choose problem-dependent
        /// parameters for the local environment. 
        /// </summary>
        /// <param name="ispec">
        /// Specifies the parameter to be returned as the value of ILAENV.
        /// <list>
        /// <item>1: the optimal blocksize; if this value is 1, an unblocked
        /// algorithm will give the best performance.</item>
        /// <item>2: the minimum block size for which the block routine
        /// should be used; if the usable block size is less than
        /// this value, an unblocked routine should be used.</item>
        /// <item>3: the crossover point (in a block routine, for N less
        /// than this value, an unblocked routine should be used).</item>
        /// <item>4: the number of shifts, used in the nonsymmetric
        /// eigenvalue routines.</item>
        /// <item>5: the minimum column dimension for blocking to be used;
        /// rectangular blocks must have dimension at least k by m,
        /// where k is given by ILAENV(2,...) and m by ILAENV(5,...).</item>
        /// <item>6: the crossover point for the SVD (when reducing an m by n
        /// matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
        /// this value, a QR factorization is used first to reduce
        /// the matrix to a triangular form.</item>
        /// <item>7: the number of processors.</item>
        /// <item>8: the crossover point for the multishift QR and QZ methods
        /// for nonsymmetric eigenvalue problems.</item>
        /// <item>9: maximum size of the subproblems at the bottom of the
        /// computation tree in the divide-and-conquer algorithm
        /// (used by xGELSD and xGESDD).</item>
        /// <item>10: ieee NaN arithmetic can be trusted not to trap.</item>
        /// <item>11: infinity arithmetic can be trusted not to trap.</item>
        /// </list></param>
        /// <param name="name">The name of the calling subroutine, in either upper case or
        /// lower case.</param>
        /// <param name="opts">The character options to the subroutine NAME, concatenated
        /// into a single character string.  For example, uplo == 'U', trans == 'T', and diag == 'N'
        /// for a triangular routine would be specified as opts = 'UTN'.</param>
        /// <param name="n1">Problem dimension for the subroutine <paramref name="name"/>. May not  be required.</param>
        /// <param name="n2">Problem dimension for the subroutine <paramref name="name"/>. May not  be required.</param>
        /// <param name="n3">Problem dimension for the subroutine <paramref name="name"/>. May not  be required.</param>
        /// <param name="n4">Problem dimension for the subroutine <paramref name="name"/>. May not  be required.</param>
        /// <returns>&gt;>=0: the value of the parameter specified by <paramref name="ispec"/>. &lt; 0:  the k-th argument had an illegal value.</returns>
        int Ilaenv(int ispec, string name, string opts, int n1, int n2, int n3, int n4);

        /// <summary>
        /// Computes a QR factorization of a real M-by-N matrix A: A = Q * R.
        /// </summary>
        /// <param name="m">The number of rows of the input matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the input matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix.
        ///  On exit, the elements on and above the diagonal of the array contain the min(<paramref name="m"/>,<paramref name="n"/>) by <paramref name="n"/> 
        ///  upper trapezoidal matrix R (R is upper triangular if <paramref name="m"/> &gt;= <paramref name="n"/>); the elements below the diagonal,
        ///  with the array <paramref name="tau"/>, represent the orthogonal matrix Q as a  product of min(<paramref name="m"/>,<paramref name="n"/>)
        ///  elementary reflectors.</param>
        /// <param name="lda">The leading dimension of the array <parmref name="a"/>.  <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="tau">Dimension min(<paramref name="m"/>,<paramref name="n"/>).
        ///  The scalar factors of the elementary reflectors.</param>
        /// <returns><list>
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// </list>
        /// </returns>
        int Dgeqrf(int m, int n, double[] a, int lda, double[] tau);

        /// <summary>
        /// Computes a QR factorization of a real M-by-N matrix A: A = Q * R.
        /// </summary>
        /// <param name="m">The number of rows of the input matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the input matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix.
        ///  On exit, the elements on and above the diagonal of the array contain the min(<paramref name="m"/>,<paramref name="n"/>) by <paramref name="n"/> 
        ///  upper trapezoidal matrix R (R is upper triangular if <paramref name="m"/> &gt;= <paramref name="n"/>); the elements below the diagonal,
        ///  with the array <paramref name="tau"/>, represent the orthogonal matrix Q as a  product of min(<paramref name="m"/>,<paramref name="n"/>)
        ///  elementary reflectors.</param>
        /// <param name="lda">The leading dimension of the array <parmref name="a"/>.  <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="tau">Dimension min(<paramref name="m"/>,<paramref name="n"/>).
        ///  The scalar factors of the elementary reflectors.</param>
        /// <param name="work">Work array with a length of at least <paramref name="n"/>*NB, where NB is the optimal blocksize.</param>
        /// <returns><list>
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// </list>
        /// </returns>
        int Dgeqrf(int m, int n, double[] a, int lda, double[] tau, double[] work);

        /// <summary>
        /// Computes a QR factorization of a real M-by-N matrix A: A = Q * R.
        /// </summary>
        /// <param name="m">The number of rows of the input matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the input matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix.
        ///  On exit, the elements on and above the diagonal of the array contain the min(<paramref name="m"/>,<paramref name="n"/>) by <paramref name="n"/> 
        ///  upper trapezoidal matrix R (R is upper triangular if <paramref name="m"/> &gt;= <paramref name="n"/>); the elements below the diagonal,
        ///  with the array <paramref name="tau"/>, represent the orthogonal matrix Q as a  product of min(<paramref name="m"/>,<paramref name="n"/>)
        ///  elementary reflectors.</param>
        /// <param name="lda">The leading dimension of the array <parmref name="a"/>.  <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="tau">Dimension min(<paramref name="m"/>,<paramref name="n"/>).
        ///  The scalar factors of the elementary reflectors.</param>
        /// <returns><list>
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// </list>
        /// </returns>
        int Zgeqrf(int m, int n, Complex[] a, int lda, Complex[] tau);

        /// <summary>
        /// Computes a QR factorization of a real M-by-N matrix A: A = Q * R.
        /// </summary>
        /// <param name="m">The number of rows of the input matrix <paramref name="a"/>. <paramref name="m"/> &gt;= 0.</param>
        /// <param name="n">The number of columns of the input matrix <paramref name="a"/>. <paramref name="n"/> &gt;= 0.</param>
        /// <param name="a">On entry, the <paramref name="m"/> by <paramref name="n"/> matrix.
        ///  On exit, the elements on and above the diagonal of the array contain the min(<paramref name="m"/>,<paramref name="n"/>) by <paramref name="n"/> 
        ///  upper trapezoidal matrix R (R is upper triangular if <paramref name="m"/> &gt;= <paramref name="n"/>); the elements below the diagonal,
        ///  with the array <paramref name="tau"/>, represent the orthogonal matrix Q as a  product of min(<paramref name="m"/>,<paramref name="n"/>)
        ///  elementary reflectors.</param>
        /// <param name="lda">The leading dimension of the array <parmref name="a"/>.  <b>lda</b> &gt;= max(1,<paramref name="m"/>).</param>
        /// <param name="tau">Dimension min(<paramref name="m"/>,<paramref name="n"/>).
        ///  The scalar factors of the elementary reflectors.</param>
        /// <param name="work">Work array with a length of at least <paramref name="n"/>*NB, where NB is the optimal blocksize.</param>
        /// <returns><list>
        /// <item>0: successful exit.</item>
        /// <item>&lt; 0: the i-th argument had an illegal value.</item>
        /// </list>
        /// </returns>
        int Zgeqrf(int m, int n, Complex[] a, int lda, Complex[] tau, Complex[] work);


    }
}
